Analytical two-center integrals over Slater geminal functions 
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We present analytical formulas for the calculation of the two-center two-electron integrals in the 
basis of Slater geminals and products of Slater orbitals. Our derivation starts with establishing a 
inhomogeneous fourth-order ordinary differential equation that is obeyed by the master integral, 
the simplest integral with inverse powers of all interparticle distances. To solve this equation it 
was necessary to introduce a new family of special functions which are defined through their series 
expansions around regular singular points of the differential equation. To increase the power of the 
interparticle distances under the sign of the integral we developed a family of open-ended recursion 
relations. A handful of special cases of the integrals is also analysed with some remarks on simpli- 
fications that occur. Additionally, we present some numerical examples of the master integral that 
validate the usefulness and correctness of the key equations derived in this paper. In particular, we 
compare our results with the calculations based on the series expansion of the exp( — 7J-12) term in 
the master integral. 

PACS numbers: 31.15.vn, 03.65.Ge, 02.30.Gp, 02.30.Hq 
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I. INTRODUCTION 

It is a well-known fact since the landmark paper of 
Kato [1] that the exact eigenfunction VP of the Schodinger 
Hamiltonian must satisfy certain conditions at the coa- 
lescence points of the particles. These are the so-called 
cusp conditions, expressed mathematically as: 

/ <9* \ 

r%(dr~J =M«?i**(r«=0) (1) 

J \ l J ' av 

where qi are the charges of the particles, fiij is the re- 
duced mass of the particles i and j , and the subscript av 
denotes the spherical average over an infinitesimal sphere 



around n 



0. The above constraint must be satisfied 



for every single pair of particles in the system. While 
the nuclear cusp condition is naturally satisfied by the 
one-electron basis constructed from the Slater orbitals, 
electronic cusp condition appears to be a far more dif- 
ficult problem. Hill [2] analysed a simple example of a 
two-electron one-center system with the basis set taken 
as the partial wave expansion: 



N 



(2) 



where Y/ m are spherical harmonics, (0i,tpi), i = 1,2, are 
the spherical angles of the vector fi , and /„ are some ra- 
dial factors. He found that the error of the energy decays 
as ~ (L + 1)~ 3 , so a rather slow convergence is obtained. 



This sad corollary can be attributed to the fact that the 
partial wave expansion has severe difficulties in fulfilling 
the electronic cusp condition. Much faster convergence 
can be expected when the basis set is extended to include 
the ri2 factor explicitly. The latter finding is a theoretical 
underpinning for a vast family of the so-called explicitly 
correlated methods. 

Explicitly correlated calculations in quantum mechan- 
ics have a long history. The first calculations of this type 
were performed on the l 1 ^ state of the helium atom by 
Hyllcraas in his classical 1929 paper [3]. The Hyllcraas 
Ansatz for the wave function of He ground state was: 
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N 



(3) 
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where s — r± + r 2 , t = r\ — r 2 , u = r\ 2 , and r, are the 
coordinates of electrons. Using a six term wave function 
of the above form with one nonlinear parameter Hyller- 
aas obtained a result with three correct significant digits 
in the ionization energy of helium [3] . The length of this 
expansion can be increased and it is a relatively easy task 
to obtain a nanohartree accuracy. Many authors tried to 
extend the form of the Hyllcraas Ansatz. For instance, 
Kinoshita [4, 5] suggested to include negative powers of 
s and u, and Schwartz [6-8] included half-integer powers 
of the latter quantities. Several researchers [9, 10] in- 
cluded logarithmic terms e.g. Log(s) in order to satisfy 
the three-particle coalescence condition of both the elec- 
trons and the nucleus. Further extension can be done 
by considering so-called "double basis set" [11-13] in 
which each combination of powers of n, r 2 , and r\ 2 is 
included twice, but with different exponential scale fac- 
tors and no logarithmic terms. Probably the most well- 
known calculations in this basis set are those of Drake et 
al. [14], where about twenty significant digits accuracy 
on the energy was reached. Of course, this idea can fur- 
ther be extended to the "triple basis set" and so forth. 
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Important from the point of view of the present paper is 
the work of Korobov [15] who obtained a 25 significant 
digits accuracy by using Slater-type geminals, i.e. the 
wave function expanded as a linear combination of the 
functions: 

4>k = e - Qfcri - /3fcr2 -T fcri2 (4) 

where a&, jk are complex parameters which were gen- 
erated quasirandomly. Recently, Nakashima and Nakat- 
suji used a method called iterative complement interac- 
tion (ICI), described in Ref. [16], and obtained forty sig- 
nificant digits accuracy which is the highest available un- 
til now. At the end of this short survey over the helium 
atom we must admit that the exponentially correlated 
Gaussian functions (ECG) were also used with success, 
see Ref. [17]. Of course, all of the methodologies men- 
tioned above can equally well be applied to the excited 
states and properties [18, 19] of the He atom, and and 
its isoelectronic series such as H~ or Li + . These systems 
were also subjects of intensive studies in the past [20-23]. 

The first explicitly correlated calculations on a molec- 
ular system, the hydrogen molecule, were carried out in 
1933 by James and Coolidge [24] with a basis set, named 
today after them (JC), of the form: 

gri&r%r&-^-K* ( 5 ) 

where and rji are elliptical coordinates. In the ad- 
vent of computers Kolos and Roothaan used this basis 
set to obtain a microhartrce accuracy in the energy cal- 
culations [25, 26]. Later on, Kolos and Wolniewicz ex- 
tended the form of the above basis set to include the 
Heitler-London function, thereby allowing to describe the 
dissociation of the molecule properly [27]. It gave rise 
to so-called Kolos- Wolniewicz (KW) basis set. The ap- 
proaches described above were subsequently extended to 
the excited states of H 2 cf. Refs. [28-30]. During the 
past decades several authors reported calculations in the 
JC [31-33] or KW [34-38] basis sets with an increasing 
accuracy. Among other approaches ICI calculations pre- 
sented by Nakatsuji et al. [39] arc worth noticing. It is 
rather astonishing that in the field of H2 ECG calcula- 
tions were proven to be very successful and even compet- 
itive with the approaches based on Slater functions [40- 
43]. Recently, Pachucki, in his tour de force paper, de- 
rived analytical equations for the integrals over the JC 
basis set [44]. This allowed to perform calculations on 
H2 with at least fifteen digits accuracy, the highest ac- 
curacy reported until now [45]. Let us end this para- 
graph by remarking that the two-electron analogues of 
H 2 , HcH+ [46-50] and He 2 + [51-53], were also studied in 
the literature. 

The lithium atom and three-electron ions are proba- 
bly the last example when Hylleraas-type basis set could 
still successfully be applied. It was possible because an- 
alytical equations for the resulting integrals [54-56] and 
useful recursion relations [57, 58] between them are all 
known. This allowed very accurate calculations, among 



which those of King [59], Yan et al. [22, 60], and Puchal- 
ski and Pachucki [61] should be mentioned. The results of 
Hyllcraas-CI and ECG calculations for the lithium atom 
are also available [20]. The accuracy of the calculations 
for the lithium atom cannot compete with that for he- 
lium. Nevertheless, the reported energy values still agree 
excellently with the best available experimental data [61]. 
The applicability of the explicitly correlated calculations 
with the Hyllcraas-like Ansatz is narrowed dramatically 
when passing to many-center and/or many-electron sys- 
tems. Since the Hylleraas-CI and ECG are the only 
methods that can be used in practice for systems such 
as beryllium atom [62-64] , the accuracy deteriorates sig- 
nificantly. Similar situation holds for other few-body sys- 
tems, H+ [49, 65, 66], H 3 [67, 68], He 2 [69], and LiH [70]. 

For many-electron systems explicitly-correlated varia- 
tional calculations are not feasible at the present. This 
is due to the high complexity in the space and pcrmu- 
tational symmetry of the wave function. However, basis 
functions including the explicit dependence on the inter- 
electronic distance r\i can be introduced into the many- 
body theory of many-electron systems. Indeed, it was re- 
alised as early as in 1966 by Byron and Joachain [71, 72] 
and later by Pan and King [82, 83], Jeziorski, Sza- 
lewicz, and collaborators [73-78] and Adamowicz and 
Sadlcj [79-81] that the pair functions appearing in the 
energy expressions of the many-body perturbation the- 
ory (MBPT), also known as the M0llcr-Plessct pertur- 
bation theory, can be expanded in terms of explicitly 
correlated functions, provided that the strong orthogo- 
nality condition is satisfied. Since the strong orthogonal- 
ity condition is difficult to meet, Szalewicz et al. [75-78] 
suggested to weaken it without loosing the mathematical 
correctness of the theory. These early explicitly corre- 
lated MBPT approaches employed the Hylleraas basis in 
the case of calculations of Byron and Joachain [71, 72] on 
the beryllium atom, and explicitly correlated Gaussian 
functions in case of the calculations on the Be, LiH, Ne, 
and H 2 systems [84-87]. In the early 1980's explicitly- 
correlated Gaussian geminals were used with success by 
Jeziorski and Szalewicz in the coupled cluster (CC) cal- 
culations [78] . One important drawback of the approach 
summarized above is that the perturbation theory and 
coupled cluster calculations involving explicitly corre- 
lated basis functions require calculations of three and 
in some cases four-electron integrals. This makes this 
kind of calculations prohibitively expensive and limited 
to small systems. A breakthrough in this respect was sug- 
gested by Klopper and Kutzelnigg [88-90] for the MBPT 
calculations and by Noga and collaborators [91, 92] for 
the CC calculations. These authors suggested to include 
only terms linear in the interelectronic distance ri 2 and 
use an approximate resolution of identity to approximate 
many-electron integrals with the two-electron integrals. 
In this way the problem of calculating many-electron in- 
tegrals was eliminated, although only in an approximate 
way. Still, this approach was shown to be very successful 
in many spectroscopic and chemical applications. See, 
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for instance, Ref. [93] for a review. Finally, the most 
recent advance in this field are the so-called explicitly- 
correlated CC-F12 methods [94-99], in which the inter- 
electronic distance, ri2, is explicitly introduced into the 
pair functions through the exponential correlation factor 
exp(— 77*12). The F12 methods have recently been im- 
plemented in an efficient manner [100-102] and shown to 
accelerate the convergence towards the basis-set limit for 
a number of properties [103-105]. Unfortunately, the F12 
method fails to reproduce accurate interaction potentials 
of diatomic molecules [106, 107], although it was shown 
to work well in the Li+LiH case [108]. 

In this paper we introduce a new basis set for accurate 
calculations on diatomic molecules, the basis of Slater 
gcminals. This basis can be used both in the variational 
calculations and in the many-body MBPT/CC theories. 
The Slater gcminal basis has several advantages over the 
explicitly correlated basis sets used in molecular calcula- 
tions thus far. Among others, it satifies both the electron- 
nuclei and electron-electron cusp conditions. Similarly as 
for atoms, the exponential correlation factor is expected 
to improve the convergence of the short-range correla- 
tions, while the Slater type one-electron part will greatly 
reduce the size of the expansion, thus leading to results 
much more accurate than possible at present. This is es- 
pecially important for the new emerging field at the bor- 
der of chemistry and physics, ultracold molecules. See the 
2012 special issue of Chemical Reviews, and in particu- 
lar papers by Quemcner and Julienne [109], Wcidemuller 
and collaborators [110], and by Koch and Shapiro [111]. 
To better appreciate the importance of high quality basis 
sets for molecular calculations on diatomic molecules, let 
us just quote calculations on the Sr2 molecule [112, 113], 
which are currently used in the interpretation of the ex- 
perimental data for the determination of the time vari- 
ation of the electron to proton mass ratio [114, 115]. 
Another very appealing application of the Slater gemi- 
nals for diatomic molecules are the calculations of the 
relativistic effects. Indeed, when the relativistic correc- 
tions are calculated in the framework of the perturbation 
method and with the Brcit-Pauli Hamiltonian it is neces- 
sary to calculate integrals with the 1/7*12 factor. Analyti- 
cal calculation of such integrals in the two-center case was 
impossible until now. It was necessary to use the infinite 
expansion in the Gegenbauer polynomials, according to 
the scheme advocated by Wolniewicz [36] . Using our ana- 
lytical equations for the integrals over the Slater geminal 
basis, all the necessary relativistic integrals involving the 
1 /r\ 2 factor are obtained by a simple one-dimensional nu- 
merical integration. Similar scheme was recently success- 
fully applied to calculation of the relativistic corrections 
for the lithium atom [61]. 

The paper is organised as follows. In Sec. II we define 
the master integral, /(r), which will serve as a generating 
integral for the calculations of all the integrals from the 
family (7) and derive a differential equation satisfied by 
/(r). In Sec. Ill we show how solve the homogeneous 
differential equation, thereby involving a new family of 



special functions. In Sec. IV we derive solutions of the 
inhomogeneous differential equation so that an analyt- 
ical expression for the master integral becomes known 
explicitly. In Sec. V we establish a family of recursion 
relations that allow calculations of the integrals with ar- 
bitrary powers of all electron-nuclear distances. Similar 
procedure is adopted in Sec. VI to let arbitrarily grow 
the power of 7*12 in the integrals. In Sec. VII we consider 
a handful of special cases of the integrals that cannot be 
calculated with the results of the previous sections. In 
these special cases, an analytical equation for the master 
integral is found in terms of well-known special functions. 
In Sec. VIII we present some numerical examples of the 
master integral that validate the usefulness and correct- 
ness of the analytical equations derived in this paper. In 
particular, we compare our results with the calculations 
based on the series expansion of the exp(— 77*12) term in 
the master integral. Finally, in Sec. IX we conclude our 
paper. 

In the paper we highly rely on the known special func- 
tions to simplify the derivation and the final formulas. 
Our convention for all special functions appearing be- 
low is the same as in Ref. [116]. We also use Mcijcr 
G-function which is defined according to Ref. [117]. 

II. THE MASTER INTEGRAL 

In this paper, we consider analytical calculation of the 
two-electron integrals in the basis of Slater geminals and 
Slater functions for a diatomic molecule. The latter basis 
set has the general form: 

0(fi , r 2 ) = r\ A r\ B r\ A r\ B r™ 2 

X e ~ U 3 r lA-U2riB-1i!2r2A-W 3 r2B-Wir 12 ' 

so it gives rise to the following class of two-electron two- 
center integrals: 

f n (i,j,k,l;u2,u 3 ,w 2 ,W3,w 1 ) = 

J d 3 7*i J d 3 r 2 r l 1A r{ B rl A r l 2B r^ 2 (7) 

X e~ UaTlA ~ u 2'"ib — ■u>2r 2A -Wsr 2 B—w l r l2 

where we adopted the following notation: r%, i = 1, 2, de- 
notes the coordinates of the electrons and r^, K = A, B, 
denotes the coordinates of the nuclei. Consequently, 
fiK = l^i — rx\ and 7*12 = \r\ — r%\ denote the electron- 
nucleus and interelectronic distances, respectively. The 
above notation will be used throughout the paper. 

It is noteworthy that the requirement 112 > 0, 773 > 0, 
t/72 > 0, 7/73 > 0, and w\ > is sufficient but much too 
strong to make the functions (6) square-integrable. This 
requirement can be significantly weakened by demanding 
only 7t2 +U3 + W1 > and 77*2 + 7773 + w\ > 0. Therefore, 
some of the nonlinear parameters can be negative without 
violation of the square- integrability principle. This result 
is reminiscent of the three-body Hylleraas integrals which 
will be discussed later. 
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If the basis set is chosen in terms of spherical harmonics 
multiplied by the radial factor and the exponential cor- 
relation factor, then using simple manipulations based 
on the ordinary trigonometric relations, one can express 
the resulting integrals in terms of combinations of the 
integrals from the family (7). 

When performing calculations for a two-electron and 
diatomic system described by the Schrodinger Hamilto- 
nian in the basis set defined by Eq. (6), all the matrix el- 
ements of the operators are readily expressed through the 
integrals (7) except for the kinetic energy operator. To 
express the latter quantities through the combinations of 
the integrals from the family (7), a somehow long deriva- 
tion is required. Not to disturb the consistency of the 
paper, this derivation is reported in the Appendix A. As 
a result, the matrix elements of the Schrodinger Hamilto- 
nian and all the integrals appearing in the nonrelativistic 
molecular physics in the basis (6) are expressed fully an- 
alytically. 



A. Definition and the momentum space 
representation 

The master integral is defined as the simplest two- 
electron integral with inverse powers of all electron- 
nuclear and intcrclectronic distances, namely: 



ri2, P2 = r 2 A, P3 = r 2 B, as: 

, d 3 p! fd 3 p 2 f d 3 p 3 e -"3Pl2 e -«2P31 



4tt J 4tt J 4ir p± 2 P31 

-W 2 p 2 p—W3P3 p —W\p\ p — Ml P23 



P2 P3 Pi P23 



(11) 



with pi2 = \p\ — p2\ and analogous formulas for p\j, and 
P23- The above representation is familiar as it is the 
generating integral from the theory of three-electron one- 
center integrals [54, 55]. Let us recall the momentum 
space representation of g(u\): 



g( Ul )=G(l, 1,1, 1,1,1), 



(12) 



where 



G(mi,m 2 ,m 3 ,m4,m 5) m 6 ) = 



1 



1 



1 



1 



1 



(13) 



B. Differential equation in the momentum space 



f(r) = 



d 3 ri f d 3 r 2 e" 



4tt 



"3 riA e — "2 r 1B 



4ir ri A 



riB 



-w 2 r 2 A P -u>3 r 2B n — wx r 12 



T2A 



T2B 



T\2 



In this subsection wc establish a differential equation 
(g) for G(l, 1, 1, 1, 1, 1). Let us first denote the integrand in 
Eq. (13) by G with an analogous notation for its param- 
eters: 



where the notation for all appearing quantities is the 
same as in Eq. (6) and r = tab is the internuclear 
distance. The reason for the choice of the multiplicative 
and the particular notation for the nonlin- 



constant 



(4tt) 



ear parameters will be clear from the further derivation. 
Once this integral is known analytically, all integrals f n 
of Eq. (7) can be obtained by multiple differentiations 
of Eq. (8) over the nonlinear parameters u 2 , M3, w 2 , W3, 
and w\. 

Our first task is to derive an analytical equation for 
the above integral. We perform a Laplace transform of 
the master integral with respect to r and therefore define 
another integral g(ui): 

g( Ul ) = J™ drf(r)e~^ r = Jfl-l^l e -^r^ (9) 

This equality allows us to calculate f(r) from the inverse 
Laplace transform formula:intcgral 



f(r) 



2ni 



ZOO + e 



dm 5(Mi)e Ul ' 



(10) 



— zoo+e 



The explicit form of the integral g(ui) can conveniently 
be written, after the simple interchange of variables p\ = 



G(mi, m 2 ,TO 3 ,m4,TO 5 ,m 6 ) 

I I I 

{k\ + uf) mi {hi + ul) m2 (fcf + uf) ms 
1 1 1 

x W2 + O mi W 3 + ^l)" 15 (Mi + ™§) m6 



(14) 



Our derivation is based on the so-called integration by 
parts identities [118, 119] and the fact that due to Green 
theorem the following family of integrals vanish: 







8^ / ^ 



d 3 kn 



d 3 fc,V., 



fc,G(l, 1,1, l,l,i; 



(15) 



where the i and j indices can independently take val- 
ues 1,2, and 3. The above identity provides nine equa- 
tions that relate the values of G with different arguments. 
These equations can be divided into three sets, the first 
set being /13, J23, -Z33 and the two other obtained by a per- 
mutation of the second index. It can be proven that to 
derive the desired differential equation only one of these 
sets has to be considered and the results from the oth- 
ers are identical. Therefore, we will consider the trio 
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I±3 , 123 , I33 but this choice is arbitrary. To give an exam- 
ple we will show the derivation for J 13 . It follows from 
the definition that: 

V 3 - [hG(l, 1,1, 1,1,1) 

- 2k! ■ k 3 G(l, 1, 2, 1, 1, 1) - 2k x ■ k 32 G{l, 1, 1, 2, 1, 1) 
+ 2^-^3(5(1,1,1,1,2,1). 

(16) 

The scalar (dot) products of several k vectors appearing 
in the above equation are expanded using the relation 

ki-k 3 = - 



h 2 - h 2 

re 13 "U 



^3 



and similar for other possible 

combinations. This allows to rewrite the r.h.s. of Eq. 
(16) as: 



k 2 - kf 
K 13 ft l 



"13 



P 



(5(1,1,2,1,1,1) 

kl + kf\ (5(1,1,1,2,1,1) 



(17) 



fc 2 3 -fc 2 -£; 2 G(l, 1,1, 1,2,1) 



The next step is to make all the coefficients multiplying 
the different G functions independent of the k vectors. 
The latter are absorbed into (5 in the following way: 

fc 2 G(l, 1, 2, 1, 1, 1) = (5(0, 1, 2, 1, 1, 1) - ulG(l, 1, 2, 1, 1, 1) 

(18) 

After necessary simplifications the expression for ii 3 be- 
comes: 

I 13 ={u 2 3 + u\-w 2 2 )G{\, 1,2, 1,1,1) 



tr, 



7/0 



W 2 )G(1,1,1,2,1,1) 



+{u 2 -u\- w 2 2 )G(l, 1, 1, 1, 2, 1) + G(l, 1, 2, 1, 0, 1) 
-G(0, 1, 2, 1, 1, 1) + G(l, 1, 1, 2, 0, 1) - G(l, 1, 1, 2, 1, 0) 
-G(l, 1,0,2, 1,1) + G(1,0, 1,2, 1,1) -G(l, 1,0, 1,2,1) 
+G(0, 1,1, 1,2,1). 

(19) 

In a very similar way the expressions for I23 and I 33 can 
be derived. The final equations are: 

7 23 =(u 2 +u 2 -w 2 )G(l, 1,2, 1,1,1) 

+{wl+ul-u\-w\)G{l, 1,1, 1,2,1) 



+(u 2 3 - wf - u 2 2 )G(l, 1, 1, 2, 1, 1) + G(l, 0, 1, 2, 1, 1) 
-G(l, 1,0,1, 2, 1) + G(0, 1,1,1,2,1) -G(l, 0,2, 1,1,0) 
-G(l, 1,0,2, 1, 1) + G(l, 1,2, 0,1,1)-G(1, 1,1, 1,2,0) 
+G(1, 1,1, 0,2,1), 

(20) 



By an inspection of these three equations we note that 
all the G integrals fall into three classes. The first class 
consists of integrals with one of the parameters mi, ttlq 
equal to zero. It is easy to verify by a direct calculation 
that these integrals belong to the class of the well-known 
Hylleraas-type helium (three-body) integrals: 



r(ni,n 2 ,n 3 ;a, /3,7) 



x e 



4- 

ari-/ 



d 3 r 2 

4-7T 
^2—7^12 



ni-1 n 2 -l n 3 - 

1 1 ' 9 '19 



(22) 



and analytical equations for these integrals are all known 
since they can be obtained from the generating integral: 



r (0,0,0; a, ^ 7) = 



1 



(0 + ^(0 + 7)08 + 7)' 



(23) 



by a proper differentiation or integration with respect 
to the nonlinear parameters a, /3, 7. Recursion rela- 
tions that enable generation of T with arbitrary values 
of m, n 2 , ri3 were presented long time ago by Kolos and 
co-workers [120]. An analytical expression to generate 
the integral T(0,0, 0) was derived earlier [121]. 

The second class of integrals consists of 
(7(1,1,2, 1,1,1), (3(1, 1,1, 2, 1,1), and (7(1,1,1,1,2,1), 
and the third class is the master integral G(l, 1, 1, 1, 1, 1). 
Therefore, we solve the set of three equations (19)— (21) 
with respect to one of the integrals from the second 
class. Let us choose G(l, 1, 1, 2, 1, 1). The result is: 



1 d 

--^G(l, 1, 1, 1, 1, 1) - 2 Wl aG(l, 1, 1, 2, 1, 1) 

2 owi (24) 



+ P(w 1 ,ui;w 2 ,u 2 ;iU3,U3) = 0, 



where a is a polynomial in all nonlinear parameters: 



22/2 2 2 1 2 
a =u 1 w 1 [u 1 — u 2 — u 3 + w 1 

„.2„„2 I „,2 1 2 „,2 



122/ 2 2,2 2 2,2\ 

+U3W3 [-U X - U 2 + U 3 - UJj - W 2 + W3) 

,222, 222, 222, 222 

+u 1 u 2 w 3 + u 1 u 3 w 2 + u 2 u 3 w 1 + w 1 w 2 w 3 . 



(25) 



J 33 =2 M 2 G(1, 1, 2, 1, 1, 1) + (w 2 -u 2 + u 2 )G(l, 1, 1, 1, 2, 1) 
+(u 2 3 +wl- u 2 2 )G(l, 1, 1, 2, 1, 1) + G(l, 0, 1, 2, 1, 1) 
-G(1,1,1,1,1,1) + G(0, 1,1, 1,2,1) 
-G(l,l, 0,1, 2,1)- (3(1, 1,0,2,1,1). 

(21) 



The function P(w\,u\\ w 2 , u 2 ; w 3 , u 3 ) is a combination of 
integrals from the first class with coefficients being some 
polynomials in the nonlinear parameters. Its derivation 
is long and does not present any advance over already 
published formulas [44, 122], so we list here only the final 
equation: 
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P(w 1 ,u 1 ;w 2 ,u 2 ;w 3 ,u 3 ) = 

= -U1W1 [(ui + w 2 ) 2 - u 2 ] T(0, 0,-1; Mi + w 2 ,u 3 , u 2 + w\) 

- u 1 w 1 [(ui + u 3 ) 2 — w 2 ] T(0, 0,-1; ui + u 3 , w 2l w\ + w 3 ) 

u 1 w 1 + w 2 u; 2 — w 3 u> 3 + ^1^2(^1 + u 2 - w 3 )] r(0, 0, —1; w\ + w 2 ,w 3 ,ui + u 2 ) 
+ [u 2 w 2 — u\w\ + + wxw 3 {u\ +u\ — wl)] T(0, 0, —1; w 1 + w 3 , w 2l ui + u 3 ) (26) 

- [^2(^2 + wi)(uf + ul - w 2 2 ) - 1*3(^1 + u 2 — w 3 )] T(0, 0,-1; u 2 + wi,u 3 ,ui + w 2 ) 

- ["3(^3 + wi)(u\ + ul - w\) — u\{u\ + ul — iu 2 )] T(0, 0,-1; u 3 + w 1} u 2 ,ui + w 3 ) 
+ wi [w 2 (u\ — u\ + wl) + w 3 {u\ + w 2 — ul)] T(0, 0, —1; w 2 + w 3 , wi 7 u 2 + u 3 ) 

+ wi [u 2 (u\ — w\+ ul) + u 3 (u\ + u\ - wl)] T(0, 0, —1; u 2 + u 3 , wi,w 2 + w 3 ), 

I 



where 



r(0,0,-l;a,/3, 7 ) 



Log 



VT+/ 3 J 



(27) 



(a-/3)(a + j a)' 

The above identity can be checked with, e.g. Ref. [123]. 
Finally, after observing that the following identity holds: 

1 dg 



G(l, 1,1, 2,1,1) 



(28) 



2wi dw\ 

one arrives at the form of the differential equation obeyed 
by g in the momentum space: 

dg 1 da . 

o"n h --s — g(ui) + P(wi,ui;w 2 ,u 2 ;w 3 ,u 3 ) = 0. 

owi 2 owi 

(29) 

By exchanging the indices at the k vectors in the defi- 
nition of G(l, 1, 1, 1, 1, 1) one can obtain analogous dif- 
ferential equations with respect to other variables. In 
particular, in the derivation the following one will be re- 
quired: 

dg 1 da 

o-z V — g(u\) + P(ui,wi- 1 u 3 ,w 3 ;w 2 ,u 2 ) = 0. 

am 2 au\ 

(30) 

The latter two equations were recently presented by 
Pachucki [44]. The solution of this differential equation 
was given by Fromm and Hill [54] and subsequently sim- 
plified considerably by Harris [56]. Unfortunately, the 
explicit form of g in terms of well-known special func- 
tions is too complicated to perform the inverse Laplace 
transform directly and obtain the two-center integrals as 
inEq. (10). Therefore, the differential equation approach 
seems to be the only way to derive analytical equations 
for the integrals family (7). 



C. Differential equation in the position space 

At this point we will depart from the previous works. 
To obtain a differential equation for the master integral 
f{r) we have to perform the inverse Laplace transform of 
the Eq. (30). Pachucki [44] performed such an inversion 



in the case of u>\ = 0, so any connection with the geminal 
basis was lost. Our case requires a generalization to an 
arbitrary physically acceptable but nonzero value of ui\ . 
Let us first rewrite the polynomial a in the following 
(convenient) way: 



a = w\u\ + fliu 2 + O2, 



o 22 22,22 22,22 

ill = ^u 2 w 1 — u 2 w 2 + u 2 w 3 — u 3 w 1 + u 3 w 2 

2 2,4 2 2 2 2 

— u 3 w 3 + w\ — w 1 w 2 — w 1 w 3 , 



n 42,222 222 222 

Si2 = u 2 w 2 + u 2 u 3 w 1 — u 2 u 3 w 2 — u 2 u 3 w 3 

222, 24 222, 42 

— U 2 W 1 W 2 + u 2 w 2 — u 2 w 2 w 3 + u 3 w 3 

222 222,24, 222 

— u 3 w 1 w 3 — u 3 w 2 w 3 + u 3 w 3 + w 1 w 2 w 3 , 



so that 



- — = Aw\u\ + 2fiiMi. 
du\ 



(31) 



(32) 



(33) 



(34) 



By inserting the above identities into Eq. (30) and col- 
lecting terms multiplying g{u\) and we get: 

(w{u{ + n lU \ + n 2 ) + (2w\u\ + Oi-ui) g(ui) 



+ P(u 1 ,w 1 ;u 3 ,w 3 ;w 2} u 2 ) = 0. 



(35) 



The inverse Laplace transform of this equation leads to: 
w\rf (i \r) + 2w{f {s) (r) + fiir/"(r) 



+ fii/'(r) + fl 2 rf(r) = U(r;w 1: u 2 ,u 3 ,w 2 ,w 3 ), 
where 

U(r;wi,u 2l u 3 , w 2 , w 3 ) = 
— du x P(u 1 ,w 1 ;u 3 ,w 3 ;w 2 ,u 2 )e UlT . 



2-ni 



(36) 



(37) 



The explicit form of U(r) is obtained by using several 
Laplace transform identities and reads: 
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with 



and 



U(r) = Y,c t U i (r) + Y^U i (r), 



i=5 



Cl = - [w 2 (u\ -u\- w\) + u 3 (wl -w% + wf)] 



C2 = - [w 3 (ul -uj + wf) + u 2 ( 



Wo — WZ — W 



I)] 



c 3 = - [w 2 (uj -ul- wl) - u 3 (wl ~wj+ wl)] , 
c -i = \ N (wj - wl- w\) - w 3 (ul - ul + wl)] , 



Ux(r) = e r ( U3 - W2 >Ei [-r (w x +u 2 + u 3 )\ - e^'^Ei [-r ( Wl +w 2 + w 3 )] 
U 2 (r) = e '-(«3-«2) Ei [_ r ( Wi +W2+ W3 )] _ e r(«2-™ 3 ) Ei [_ r ( Wi +u . 2 + ug)] 

U 3 (r) = e -K«3+»2) { Ei [_ r (wi + U2 _ Ua j] _ Ei [_ r ( U2 _ U3 _ W2 + W3 j] 
+ Ei [-r (u>i - w 2 + u> 3 )]} - e^^+^^Ei [-r («a + u 3 + w 2 + w 3 )] 



+ e - r{u3+W2) Los 



(wi +w 2 + w 3 ) (u 2 + U 3 + Wl) (u 2 - u 3 - w 2 + w 3 ) 



(w\ - w 2 + w 3 ) (wi +u 2 - u 3 ) (u 2 + u 3 + w 2 + w 3 ) 



U 4 {r) =e~ r{u2+W3) {Ei [-r (w x + u 3 - u 2 )] - Ei [-r (u 3 - u 2 - w 3 + w 2 )} 
+ Ei [-r (tui - w 3 + tu 2 )]} - e r (" 2+,i ' 3) Ei [-r (it 2 +u 3 + »2 + to 3 )] 



+e -r(u 2 +«, 3 ) Lo 



(ii'l + W 2 + W3) («2 + "3 + Wl) (u 2 - U 3 — W 2 + W 3 ) 



(wi — w 2 + w 3 ) (wi + u 2 - u 3 ) (u 2 + u 3 + w 2 + w 3 ) 



(38) 

(39) 
(40) 
(41) 
(42) 

(43) 
(44) 

(45) 



(46) 



U 5 (r) = - !^i e -K«2+»3) (I + U2 + w \ + H!l e -r(u 3 + Wl + W3 ) (l +Ua + Wl + W3 
r \r I r \r 



_ e -r(u 2 +w 3 ) Wi [ 7 + + (1 _ 7 ) £'( r )] + w 2 e -r(u 3 +w 1+ w 3 ) 

+e -r( U3 + Wl+W3 ) [7 (wg + W3+ Wl)S (r) + (1 - 7) 6'(r)] , 



- + 7<5(r) 
r 



(47) 



U, 



;(r) = - ^ e -K«3+™ 2 ) ( I + U3 + W2 ) + HH e -r(u 2+Wl+W2 ) (l + U2 + Wi+w \ 
r \r J r \r J 

_ e -r(u 2 +w 3 ) Wi [7 (u3 + W2) ^ + (i _ 7 ) s'(r)] + w \e- r ^ +w ^ +w ^ 
+e -r(u 2+Wl +w 2 ) [7 (W2 + W2 + Wi) ^ + (1 _ 7 ) S '(r)] , 



- + 7<5(r) 
r 



(48) 



C/ 7 (r) £(r) Log(w;i +w 2 + w 3 ), U s (r) = w\ S(r) Log(wi +u 2 + u 3 ). 



(49) 



r 



The results presented above require some comments. 
First of all let us establish the connection with the 



Pachucki differential equation, the zero limit in w\ of 
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Eq. (36). By setting w\ = and observing that: 

2__^2, jp/ \ _ U(r;0,u 2 ,u 3 ,w 2 ,w 3 ) 

P ~ fa Ll=0 ' { ' ~ fii 



one arrives at: 

rf"(r) + f\r) - p 2 rf(r) + F(r) = 0, 



(50) 



(51) 



which exactly coincides with the result given by 
Pachucki [44] for the case of orbital basis. Second, at 
the end of this subsection we would like to mention that 
in the further derivation we will make use of two other 
functions which arc obtained as the inverse Laplace trans- 
forms of P, namely: 



W(r;w 1 ,U2,U3,W2,W3) = 

— / dm P(wi,ui; w 2 ,u 2 ; w 3 ,u 3 )e Uir , 
Ziyi l /_j oc+e 



V(r;w 1 ,u 2 ,u 3 ,w 2 ,w 3 ) 

/>ioo+e 

27ri J_ ioo+e 



dm P(w 3 ,u 3 \w 2 ,u 2 ;wi,u 1 )e Ul 



(52) 



(53) 



Since explicit formulas for these functions have not been 
presented in the literature thus far, we list them in the 
Appendix B. 



III. SOLUTION OF THE HOMOGENEOUS 
DIFFERENTIAL EQUATION 

First, we will solve the homogeneous version of the 
^eminal differential equation: 



wlrf^{r) +2w 2 1 f^(r) + n 1 rf^(r)+Q 1 f' H (r) 
+ n 2 rf H {r) = 0, 



(54) 



where the subscript H was added to designate the so- 
lution of the homogeneous equation. The above equa- 
tion is a homogeneous linear ordinary differential equa- 
tion (ODE) with non-constant coefficients. We found 
it very difficult, if not impossible, to express the solu- 
tion in terms of well-known special or analytical func- 
tions. Any manipulations performed with Eq. (54) were 
proven fruitless in bringing this equation into a charac- 
teristic form, thus enabling an analytical solution. It was 
also impossible to find the solution by using a symbolic 
mathematical package such as Mathematica [124]. 

It is interesting from the mathematical point of view 
that Eq. (54) can be brought to the form 



dr 2 



dr 2 



dr \ dr 



= -rQ 2 f H (r), (55) 



so that it can be considered as a generalization of the 
Sturm-Liouville (S-L) equation to the fourth order ODE 



with the weight (or density) function equal to r and eigen- 
value — fl 2 . 

Because of all the above, we decided to define a new 
family of special functions which, by definition, form 
the general solution of the differetial equation (54) . We 
will find its form by using the generalized version of the 
Frobenius method (see, e.g. Rcf. [125]). Precisely, we will 
find a solution in terms of the series expansion around 
two singular points, zero and infinity. Our first Ansatz is 
an ordinary regular expansion around r = 0: 



atr 



(56) 



fe=0 



We insert this Ansatz into Eq. (54), collect terms multi- 
plying the same power of r and require them to zero to 
make the differential equation satisfied for all values of 
r. This establishes the recurrence relation that connects 
the values of a k with different k. The final result reads: 



w\{k 
tt 2 a,k- 



l)(k + 2) 2 (k + 3H+3 + fii(fc + l) 2 a k+1 + 
= for k > 1, 



and the indicial equation: 

12w\a 3 + f^idi = 



(57) 



(58) 



Equations (57) and (58) need to be simultaneously sat- 
isfied. However, there is a freedom in the choice of three 
initial parameters a ,ai and a 2 . Therefore, we specify 
three new special functions Li(r), i = 1,2,3, using their 
expansions around r = given by Eq. (56) and the re- 
currence relation (57). The choice of the three initial 
parameters is conventional and we put: 



Li(r) with clq = 1, oi = 0,a 2 = 0, 
L 2 (r) with do = 0, a± = 1, a 2 = 0, 
L 3 (r) with ciq = 0,ai = 0,a 2 = 1. 



(59) 



This convention will be used throughout the paper. Let 
us justify the choice of the formulas (59). One may argue 
that choice a = 1 in L\ is very special but by putting 
so = C ^ 1 we obtain a function which is just L\ mul- 
tiplied by C. The choice of a multiplicative constant is 
immaterial in our context and, consequently, so is the 
choice of C . The same is true for the values of a\ and a 2 
in L 2 and L 3 , respectively. Similarly, by defining a func- 
tion with oo = ljOi = 1, for example, we obtain a linear 
combination of L\ and L 2 . Because of these properties, 
we find the convention (59) justified. 

It is clear that the three functions obtained in the pre- 
vious paragraph are not sufficient to give the general so- 
lution of the homogeneous geminal differential equation. 
Our second trial for the expansion around r = is some- 
what less obvious: 



f H {r)=L l {r)hog(r) 



k=0 



b k r k 



(60) 
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where the coefficients b k are to be determined by insert- 
ing the expression (60) into the homogeneous differen- 
tial equation and collecting terms multiplying r k and 
r fc Log(r). This results in the recurrence relation: 

2n iai + 28wja 3 + fii&i + 12wjb 3 = 0, (61) 
2kfl iak + 2{k + l)(2k 2 + 4k + l)wla k+2 + fl 2 b k _ 2 
+ fc 2 Qib fc + w\k{k + l) 2 (fc + 2)b k+2 = 0. 

(62) 

Additionally, as soon as w\ ^ the above Ansatz re- 
quires ao = 0, ai = 1, a 2 = 0. As before, we have three 
parameters which can be chosen freely, bo, b\ and b 2 . 
Since we seek for only one function let us put 60 = 0, 
b\ = 1, b 2 = which leads to: 



and for k > 2: 



U(r) = L 2 {r)Log(r) + V b k r 



k=l 



(63) 



One can show that any function constructed with a dif- 
ferent choice of bo, b\ and b 2 can be expressed as a linear 
combination of Li(r),L 2 {r),L 3 {r),L^{r). This formally 
completes the solution of the homogeneous differential 
equation (36). 

The expansions around r = presented above are con- 
vergent for all finite values of r since the coefficients mul- 
tiplying the powers of r decay faster than any polynomial 
when k — > 00. However, the rate of convergence of these 
series can be expected to be poor for large values of r and 
therefore prohibit an accurate calculation in this regime. 
As a result, it might be beneficial to obtain their asymp- 
totic expansion which will be valid and rapidly conver- 
gent for large values of r. The latter expansion can be 
constructed from the Ansatz: 



f H (r) = e tr qfcr" 



k — p 



(64) 



k=0 



where a k , t and p are coefficients to be determined. By 
inserting this trial function into the differential equation 
and grouping coefficients multiplying the same powers of 
1/r one obtains indicial equations specifying p and t: 



wjt 4 + n^ 2 + n 2 = o, 
1 



(65) 
(66) 



and the recursion relation for a k with the value of p al- 
ready fixed at 1/2: 



1 3 

-f2ia + -t 2 w\ao - 2Qiiai - U 3 w\ai = 0, (67) 
9 27 

— 3tw\ao + — Oi«i + — t 2 w 2 cii — AQita 2 — 8t 3 w 2 a 2 = 0, 

(68) 



=— wja k (2k + 3) 2 (2k + l)(2fc + 5) 
16 



-tw\a k ^ 
4(2fcH 



i(fc + 2)(2fc + 3)(2fc + 5) 
5) 2 a fc+2 (fii + m 2 w\) 



(69) 



-2(fc + 3)a fe+3 (fix* + 2t 3 wl] 



Eq. (65) has four solutions U, i = 1, ...,4, which corre- 
spond to four functions determining the general solution 
of the homogeneous differential equation. We see that 
it is dependent on the sign of ij whether convergent or 
divergent expansion is obtained. The final result can be 
written as: 



Jr f— ' r K 

v fc=0 



(70) 



where ao can freely be chosen. 

It is interesting to establish a connection between the 
new special functions Li and the modified Bessel func- 
tions of the first, Io( r ), and the second, Ko(r), kind. By 
setting wi = Eqs. (57) and (58) become: 

ni(fc + i) 2 o fc+ i + n a o fc _i = o, 01 = o, (71) 



so that the recursion can be solved explicitly to give: 

k 

u>i=Q P _ n r*7o\ 

a2k ~ 2 2k kl 2 ~ 2 2 ~ k ~k\ 2 ~ ' ~ ' ( ' 

and the series can be brought into the closed form: 



Ik 

Hi 



E 

fe=0 



^2 A; ^,2 A; 



Io(pr), 



(73) 



coinciding with the Bessel function of the first kind. Sim- 
ilarly, by setting w\ = in Eqs. (61) and (62) one finds a 
linear combination of Io(j>r) and Ko{pr) to be the Wi = 
limit of I/4(r). One could force the exact relationship: 



lim L 4 (r) = K {pr), 

uii — >0 



(74) 



by a proper choice of the initial parameters. Our choice 
was made for the sake of simplicity as indicated before. 
Similar result is found with the asymptotic expansions of 
Li(r). Whenever w\ = 0, Eq. (65) has two solutions: 



-I 



(75) 



so that Eq. (70) becomes the asymptotic expansion of Jo 
(with t — p) or Ko (with t = —p) . 

We believe that because of the interesting properties 
of the Li(r) functions and their strong connection with 
the Bessel functions they can be understood as a gener- 
alization to the fourth order differential equation. There- 
fore, we give them the name hyper-Bessel functions. In 
analogy, L\,L 2 ,L 3 functions are hyper-Bessel functions 
of the first kind and L4 is the hyper-Bessel function of 
the second kind. 
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IV. SOLUTION OF THE INHOMOGENEOUS 
DIFFERENTIAL EQUATION 

The next step in our derivation is to use the properties 
of the functions introduced in Sec. Ill to obtain solution 
of the inhomogeneous differential equation (36). In this 
work we decided to use the method based on the Wron- 
skian determinants. Starting with the general solution of 
the homogeneous equation: 

f H (r) = CiLi(r) + c 2 L 2 (r) + c 3 L a {r) + c 4 L 4 (r), (76) 

we make coefficients Ci explicit functions of r, yi(r), and 
require the combination 



f(r) = y 1 {r)L 1 (r) + y 2 (r)L 2 (r) 
+ y3(r)L 3 (r) + y 4 (r)L 4 (r), 



(77) 



to satisfy the inhomogeneous gcminal equation (36). Dif- 
ferentiation of the above equation leads to: 



f'(r)=[y 1 (r)L' 1 (r) + y 2 (r)L' 2 (r) 
+ y 3 (r)L' 3 (r)+y 4 (r)L' 4 (r)} 
+ [y[(r)L 1 (r)+y' 2 (r)L 2 (r) 
+ y' 3 (r)L 3 (r)+y' 4 (r)L 4 (r)}. 



(78) 



Requirement that the second expression vanishes identi- 
cally for all r: 



= |/i(r)ii(r)+i^ a (r)L 2 (r) 
+ y 3 (r)L 3 (r) + y' A (r)L 4 (r), 

enables us bring the first derivative into the form: 

f(r) = yi (r)L[(r)+y 2 (r)L' 2 (r) 
+ y3(r)L' 3 {r) +y 4 (r)L' 4 (r). 

Similarly, higher-order derivatives are found to be: 

f"(r)=y 1 (r)L'{(r) + y 2 (r)L%(r) 
+ y 3 (r)L' 3 ^r)+y 4 (r)L'l(r), 

f {3 \r)=y 1 {r)Lf\r)+y 2 {r)L 2 3) {r) 
+ y 3 (r)L 3 3 \r)+y 4 (r)L (3) (r), 

f (A) (r)=[y 1 (r)L ( t\r)+y 2 (r)L^(r) 
+y 3 (r)L 2 i \r)+y 4 (r)L 4 i \r) 
+ [y' 1 (r)L[ 3 \r)+y' 2 (r)L 2 3) (r) 
+ y' 3 (r)L^(r)+yi(r)Lf\r) 



(79) 



(80) 



(81) 



(82) 



(83) 



where additional constraints on yi(r) where imposed: 

y[(r)L[(r)+y' 2 (r)L' 2 (r)+y' 3 (r)L' 3 (r)+y' 4 (r)L' 4 (r) = 0, 

(84) 

y[(r)L'{(r) + y' 2 (r)L'*(r) + y' 3 (r)L' 3 \r) + y' 4 {r)L'l{r) = 0. 

(85) 



By inserting Eqs. (81), (82), and (83) into the differential 
equation and noting that the functions Li(r) satisfy the 
homogeneous differential equation one arrives at: 

U(r) =wfr [y[(r)L[ 3 \r) + y' 2 (r)L 2 3) (r) (86) 

+ y' 3 (r)L 3 3) (r)+y' 4 (r)L{ 3) (r)] . (87) 

The above equation together with Eqs. (79), (84), and 
(85) specify the four-dimensional system of linear equa- 
tions: 

4 3) (r) L {3 \r) L {3 \r) hf{r) 
L'{(r) L»{r) L»{r) L'{{r) 
L[(r) L' 2 (r) L' 3 (r) L' 4 (r) 
L x (r) L 2 {r) L 3 (r) L 4 (r) 



rvi(r)- 




- U{r) 
w\r 


y' 2 (r) 







y' 3 (r) 







ly'Ar). 








that can easily be solved symbolically for y[{r) by using 
the Cramer's rule. The result reads: 



y[(r) 



U(r) W 4 {r) 
w\r W{r) ' 



(89) 



where W{r) is the Wronskian determinant of the func- 
tions Lx(r), L 2 (r), L 3 (r), L 4 (r): 



W{r) = 



L»{r) 
L'i{r) 
Li(r) 



4 3) ( 

L'i{r) 
L' 2 {r) 
L 2 (r) 



-(3) 



L' 3 \r) 
L' 3 {r) 
L 3 (r) 



, (r) I*™ 

L'l(r) 
L' 4 (r) 
L 4 (r) 



(90) 



and Wk{r) are the same as W(r) apart form the fact that 
thefcth column was replaced by the unit vector [1, 0, 0, 0], 
for example: 



W 4 (r) 



1 L 2 3) (r) Lf(r) Lf(r) 

LZ(r) L»(r) L'l(r) 

L' 2 (r) L' 3 (r) L' 4 (r) 

L 2 (r) L 3 (r) L 4 (r) (91) 

L' 2 \r) L>'{r) L'({r) ' 
L' 2 (r) L' 3 (r) L' 4 (r) 
L 2 (r) L 3 {r) L 4 {r) 



Some degree of suspicion may be connected with the 
fact that Wronskian appears in the denominator. How- 
ever, since the functions L 4 (r), L 2 (r), L 3 (r), and L 4 (r) 
span the space of solutions of the homogeneous differ- 
ential equation they cannot be linearly dependent and 
thus W(r) cannot vanish identically. Eq. (89) can now 
formally be integrated 



yi( r ) = dr 



U{r) Wx{r) 
w\r W(r) 



(92) 



so that the solution of the inhomogeneous differential 
equation is 



f(r) 



4 

E 

i=l 



Li{r) / dr 



U(r) Wj(r) 
w\r W(r) ' 



(93) 
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where the initial conditions have not been imposed yet. 
At this point we observe that the solution is rather com- 
plicated because of the presence of five determinants, in- 
cluding the Wronskian itself which is the most cumber- 
some in the calculations. Therefore, it will be advan- 
tageous to introduce some simplifications in the above 
formula. It turns out that the appearance of the Wron- 
skian can be eliminated altogether by using the so-called 
Abel's identity which states that for any nth-order ho- 
mogeneous differential equation of the form: 



f'(i) +Pn-i (x)y (n - 1) (x) + . . .+p (x)y(x) = 0, (94) 

the Wronskian constructed from n linearly independent 
solutions can be expressed as: 



W(x) = W(x )exp 



dx'p n -i(x' 



(95) 



provided that p n -\{x) is continuous on the interval 
[xo,x]. In the case of the geminal differential equation 
p n _x(r) = - (continuous on < r < oo), so that inte- 
gration can easily be carried out and the Abel's identity 
is: 



W{r) = W(r ) ( ^ 



(96) 



We need to specify the point tq. Our choice r*o = 1 was 
motivated by the fact that the above equation takes a 
very simple form and that ro = 1 is sufficiently close to 
r = at which series expansions of L-s(r) were provided. 
It allows a robust calculation of W(l) for any values of 
the nonlinear parameters. By inserting the identity: 



W(r) 



(97) 



into the solution (93) considerable simplifications occur: 



/(*•) = E 



w\W{l) 



drrU(r)Wi(r). 



(98) 



Despite a considerable effort we did not manage to sim- 
plify this equation further, at least in the general case. 
Such a simplification will occur for a special case consid- 
ered in the next subsection. 

Finally, we have to impose four initial conditions on the 
above solution to make it consistent with the definition 
of the master integral. The first three initial conditions 
are natural: 



/(0) = 0, lim f(r) = 0, lim f'(r) = 0. 



(99) 



The fourth initial condition is somehow more compli- 
cated, but we see that whenever r — > then t\b T\a — 
n etc., so that in the r = limit the derivative of the 
master integral becomes the well-known integral: 



no) 



<t'> 



4tt 



r 2 e 



-(v.3+u 2 )r 1 p -(w 2 +w 3 )r2 „-i«iri2 



4tt 



r{ r$ ri2 

(100) 

Analytical formula for this integral is well-known, cf. Eq. 
(34) of Rcf. [123]. After the four initial conditions are 
imposed the master integral becomes: 



fir) = J M-J r dr'r'U{r')W 1 {r') ^ 



w(W(l) 



w(W(l) J r 



dr' r'U(r')W 3 {r') + -^^ 



wfW(l) Jo 



dr' r'U{r')W 2 {r') 
dr' r'U{r')W 4 (r'). 



(101) 



One can check that the above formula satisfies the initial tained recursively by differentiation of Eq. (36): 
conditions (99) and (100), and therefore is the solution of 

the inhomogeneous differential geminal equation. In the f(n+4)/ \_ U^(r) n + 2 ( n+3 ), . Oi („ +2 ). 

further derivation we will also need the values of /'(r), ^ w — rw 2 r w\ 

f"(r) and P 3 **(r). Higher-order derivatives can be ob- ( m _i_ -\\ o o „ 

J W 6 + ]) f{n+l)/ r \ S '2 j( n ) ^gn ^(n-l) 

rwf ' w\ rw\ 

(102) 
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The first derivative of the master integral with respect to r is obtained directly from the representation (101): 

m - 4w> i *' r ' ' Ar ' )wi<r ' } - r * r ' u ^ w ^ 

- f *' r ' u(r ' )w ' {r,) + $m L *' *•' D(rW) < 103 > 

i 



The non-integral term is equal to: 

i— 1 * ^ ' i=X 



so that values of the latter three quantities can be cal- 
culated with an insignificant additional cost once the nu- 
merical integration of the integrals appearing in f(r) is 
done. 



V. RECURSION RELATIONS FOR THE 
POWERS OF r 1A ,r 1B ,r2A,r 2 B 

With the value of the master integral at hand, we turn 
to the calculation of the integrals with arbitrary powers 
of r\A , t\b j j f2B ■ They are obtained by differentia- 
tion of the master integral with respect to the nonlin- 
ear parameters. Explicit differentiation of Eq. (104) is 
cumbersome and connected with painful and expensive 
numerical integrations. Therefore, to start the recursion 



so it vanishes identically on the basis of the initial as- 
sumption (79). The first derivative of f(r) becomes: 



(105) 



(106) 



(107) 



I 

relations, we must establish an equation that connects 
the value of the derivative of the master integral with 
respect to, say, W3, to the master integral and optionally 
its derivatives with respect to r. The latter quantities 
can be computed by using the theory presented in the 
previous section. 

The desired recursion relation can be derived from two 
differential equations in the momentum space. The first 
was already derived in the subsection II B: 

dg 1 da 

<y-R 1- r a — g(ui) + P(ui,wi;u 3 ,w 3 ;w 2 ,U2) = 0, 

am 2 am 

(108) 

and the second is obtained by the proper exchange of the 
nonlinear parameters, making use of the fact that both 



- i§h f *' *•' u{T ' )mr ' ] + 3m L *' '' u(rWy 

Similarly, using the conditions (84) and (85), explicit formulas for f"(r) and f^ 3 '(r) are obtained: 

/ " w - 4w> f. ""' u(T ' m <T,) -^§mf *' r ' u(T ' m<T,) 

- iM) f *' r ' u(r ' )W3(r ' } + 3m L dr ' u{r ' )mr% 

'"'^ -^mi" r ' u{r ' W) - 3m I *' r ' V{r ' )w ' ir,) 
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q and a are invariant under the latter operations: 

3/(4) Qf" df ldCli , 

w i~5 — + fi 2^ — + ^-5 — f [r) 

dg i da dw * dw * dw * 2 dw * (111) 

^ + 2Q^9(ui) + P(w 3 ,u m ur,m,u 1 ) = 0. +i^/(r) + V(r;w 1 ,« a> t« 8| w a> w 8 )=0 1 

(109) 2 d ™ 3 

By taking the inverse Laplace transform of both equa- 
tions one obtains 

-v$rfW(r) -2 W 2fW(r) -Q ir f"(r) - fii/'W , x . . , , n 

(110) and by differentiation of the first equation with respect 

-fiar/(r) + E/(r; Wl , u 2 , u 3 , w 2 , w 3 ) = 0, to W3 one obtains a pair; 



Ei 



2w( 



3/(3) 



9/" 90i 



5/' 



_ r/(r ) - 



9C/(r) 



9l«3 



0, 



(112) 



a/ (4 



<9u>3 



n 5 



<9u>3 



r 



where the notation for the nonlinear parameters in U 
and V was suppressed for brevity. These two equations 
provide a starting point to establish an explicit recursion 
relation. However, its derivation is still a nontrivial task 
since Ei, E 2 , apart from the desired term consist of 
the derivatives of the latter with respect to r up to the 
fourth order. Our approach was based on the following 
three additional identities that are defined as: 

Q 

E 3 = ^(E!+rE 2 ), (113) 

E 4 = ^(rE 3 -2E 1 ), (114) 

E 5 = ^(E 4 -4E 2 ). (115) 

The reason for making the combinations above is as fol- 
lows. At each step we cancel out the fourth-order deriva- 
tive of J^- with respect to r and then create it back by 
doing a differentiation with respect to r. By repeating 
this procedure three times we figure out that the Eq. 
(113) is a set of equations with five unknown quantities: 

df df df" dfW df^ 

UW3 OW3 OW3 OW3 OW3 

so it can be solved analytically. The differentiation per- 
formed at each step guarantees that E,;, i = 1,5, are 
linearly independent as long as none of the coefficients 
multiplying the unknown quantities in the initial equa- 
tions for Ei and E 2 vanishes. Higher-order derivatives of 



0^- over r do not appear. The final result is: 

df_ = 

dw 3 

2wl^ Uv(r) + 2rV'(r) + 2^- + f^[/(r) + rf'(r)} 
D { dw 3 dw3 

" S^ [/ " (r) + r/(3) ^] j + Wo { ~ W[r) ~ 2rV ' {r) 

du>3 dW3 dW3 J 

+ W l^{-lW"(r)-2rV^{r)-2 dU[ * ){r) 
D Q { dw 3 

+ fV'M + rf^(r)} + |%/ (4) M + rf^(r)}\ , 
dw3 dw3 J 

(116) 

where Do is the common denominator: 

D = 2fl 2 (Q\ - 4wln 2 ) ■ (117) 

It is noteworthy that the procedure in which the required 
set E^ i = 1,5, was obtained is somehow ambiguous. 
Only the first step of this procedure, formation of E3, is 
unique since there is only one correct method to obtain a 
useful equation by cancelling out the fourth-order deriva- 
tive of a^— . In the further steps such an elimination can 
be performed using different equations which were ob- 
tained previously and the number of possibilities grows 
with the number of steps taken. We cannot prove that 
the particular choice of equations for Ej, i = 1,5, which 
we used here is 'the best'. However, in our procedure we 
tried to minimize the order of the derivatives of functions 
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U(r) and V(r) that appeared in the final result. It leads 
to equation (116) which turned out to regular. 

By multiplying both sides of the relation (116) by Dq 
and by further differentiation one can calculate arbitrary 
derivative over the nonlinear parameters thus advancing 
the powers of Tia,t\b,T2A-,T2b as much as necessary. 
The recursion relations for the derivatives over w 2 , "3 
and u 2 that cannot directly be calculated from the above 
formula are obtained with the use of the symmetry of 
the master integral. Namely, by permuting W2 <-t W3 and 
u 2 <-> U3 (exchange of the nuclei A B in the master 
integral) and noting that the master integral is invariant 
with respect to this permutation, analogous recursion re- 
lation for is obtained. Similarly, the exchange of 
u 2 <-> W3 and 7*3 •*-> W2 (change of the electrons' num- 
bering 1 f>2) results in the derivative over u 2 . Finally, 
the use of both of these permutations gives the derivative 
over 7/3. 

We listed only the formula for -S^- despite the fact 
that by solving the set of equations for Ej, i = 1,5, its 
derivatives over r up to the fourth-order are obtained 
as by-products. From the mathematical point of view 
equivalent formulas can be derived by differentiating Eq. 
(116) over r. Although numerical results obtained in this 
manner are the same, formulas for higher-order deriva- 
tives over r calculated from the solution of Eqs. (113) 
are much more transparent. In particular, they do not 
include higher-order derivatives of the functions U(r), 
V(r), and of the master integral. Therefore, we list all 
the missing formulas in the Appendix C. 



VI. RECURSION RELATIONS FOR THE 
POWERS OF ri2 



Since the integrals in the Slater geminal basis consid- 
ered here already include explicit correlation factor, there 
is a little point in growing powers of 7*12 in the initial ba- 
sis set. However, such a possibility is open and we will 
elaborate it in this section. In particular, we will derive 
an analytical equation for the overlap integral over Slater 
gcminals which, despite its simplicity at first glance, has 
not found analytical solution yet. Our approach is sim- 
ilar to the one presented in the previous section. We 
shall establish a relation that connects the value of -S^- 
with the master integral and its derivatives over r. In 
the derivation we use the following differential equations 
for g in the momentum space: 



dg 1 da 

o'a h o"5 — 9( u i) + P{ui,wi;u 3 ,W3;w2,U2) = 0, 



dg 1 da . . 

^77" + 7:^ 9{Ul) +P[Wl,Ui\ W 2 ,U 2 ; W3,U 3 ) = 



(118) 

: 0. 

(119) 



The first of these equations is differentiated with respect 
to wi and then the inverse Laplace transform is per- 
formed to give: 



J 



Ei 



2wf 



2 df^ 



1 dw\ 1 dw\ 
2rw 1 f i4) (r)-4w 1 f (3) (r) 



Of" 



art! 



Oi 



df 
dwi 



n 9 r 



dwi 



f'(r) 



on 2 



dU(r) 



(120) 



0. 



4) 



df 



1 dflo 



+ ^ok + wi/(4) w + 2^/'^ + 2^ /w + w{r) 



(121) 



Using a similar procedure as for the derivatives over 7JJ3 
we form a set of equations: 



which are then solved for the following quantities: 



E 3 = ^ (Ei + rE 2 ) , 
E 4 = ^ (E 3 + 2E 2 ) , 



E, 



dr 



(122) 
(123) 
(124) 



df df df" df^ a/(4) 
duii ' dwi ' dwi ' dwi ' dw\ 
The final equation for J^j- is given by: 
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df_ = 
dwi 

~2w 2 ^- i-6W(r) - 2rW'(r) - 2^- - g[/(r) - r/'(r)] 
+ +rf (3) (r)} +6 Wl f^(r) + 2w 1 f^(r) >l j 

+ ^1 f _ 4W (r) - 2rW\r) - 2^- + |^r/'(r) + |%/"(r) + r/< 3 >(r)] + (125) 
L>o L OtWl OWl OWi 

+ 8i0i/ (4) (r) +2ru;i/ (5) (r)} 

(-10W"(r) - 2rW«(r) - 2^^ + fV'M + r/( 3 >(r)] 
+ l^[3/ (4) (r) + r/^(r)] + 10wi/^(r) + 2 Wl r/ (7) (r)} . 



Higher powers of r^i are obtained by further differenti- 
ation of the above equation. As before, useful formulas 
resulting from the solution of the set for E,, i = 1,5, are 
listed in the Appendix C. 

VII. SPECIAL CASES 

In this section we consider four special cases of the in- 
tegrals corresponding to situations when coefficients fii 
and/or f2 2 vanish. From the mathematical point of view 
the recursion relations for the coefficients a k and bk re- 
main valid since in their recursive evaluation one never 
divides by fii or f2 2 . Therefore, the representation of the 
master integral given by Eq. (100) is still correct. How- 
ever, for practical reasons it is useful to consider these two 
special cases in details since the solution of the homoge- 
neous differential equation can be expressed in terms of 
well-known special functions. This makes the implemen- 
tation of the method much simpler. 

A. Vanishing Q2 coefficient 

Vanishing H 2 coefficient is probably the most impor- 
tant special case since it occurs for a handful of physi- 
cally important classes of integrals. This includes expo- 
nentially correlated analogue of the symmetric James- 
Coolidge basis set [24] (112 = 113 = W2 = W3 = x) 
and symmetric exchange integrals over atomic orbitals 
(it3 = w>2 — x, U2 = W3 = y). Singularities also appear 
whenever: 

w\ = {u 2 2 wl - u\wl) ( 2 1 2 - 2 1 g ) . (126) 

Moreover, the recursion relations established in the previ- 
ous subsections are not valid in this case since f2 2 appears 
in the denominator in the key formulas. 



In this special case the geminal differential equation 
takes the form: 

w\rfW{r) + 2wl.f^{r) + n ir f"(r) + n x /'(r) = U(r), 

(127) 

so that the homogeneous equation is: 

wlrffir) + 2wlf { *\r) + fiir/£(r) + faf^r) = 0. 

(128) 

The simplest way to obtain the solution of the latter 
equation is to use the recursion relations for the coef- 
ficients in the series expansions that were derived for the 
general case, Eqs. (57), (58), (61), and (62). By setting 
f2 2 = they become: 

w 2 1 (k + 2) 2 {k + 3)a k+ 3 + n 1 (k + l)a k+1 =0 for k > 1, 

(129) 

12wja 3 + fiioi = 0. (130) 

As before, there are three initial parameters that can 
freely be chosen, ao, 01, a 2 . Let us make the same choice 
as in Eq. (59): 

L\{r) with ao = l,oi = 0,a 2 = 0, 

L 2 {r) with ao = 0, a\ — 1, a 2 = 0, (131) 

L>3(r) with ao = 0, 01 = 0,a 2 = 1. 

The solutions of Eq. (128) are denoted by tilde 
to distinguish them from the solutions in the gen- 
eral case. The resulting functions can be expressed 
in terms of the generalized hypcrgcometric function 
(pF q [on, ... , a p ; 61, ... , b q \ z\) and some elementary func- 
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tions: 



li(r) = 1, 
L 2 (r) = r iF 2 

L 3 (r) = r 2 2 F 3 



1 3 

2' ' 2' 



' '2'2' ' W ? 



(132) 
(133) 

(134) 



The recursion relation for the coefficients b k becomes: 

2fl iai + 2%w\a 3 + ilib! + 12wlb 3 = 0, (135) 
2ktt 1 a k + 2(k + l)(2k 2 + Ak + l)w\a k+2 

(136) 

+ k 2 n x b k + w\k{k + i) 2 (k + 2)b k+2 = o. 

At this point it is very useful to depart slightly from the 
previous approach and choose a little less obvious initial 
conditions for the scries b k : 



b Q = 0, h = - 
ir 



Log 



2wi 



7-1 



b 2 = 0. 



(137) 



With this choice the function L4 (r) takes a very appeal- 
ing form: 



U{r) 



Y 

Wi 



il'i 



r)H_ 1 



Wl 



r Ho 



/H7 



Wl 



(138) 
(139) 



where Y a is the Bessel function of the second kind and 
H a is the Struve function, both of the order a. This 
completes the solution of Eq. (128). In this particular 
case we found a closed expression for /s(r) in terms of 
the known special functions, so that the implementation 
and numerical realisation becomes significantly simpler. 
Since the initial conditions for /(r) in this special case 
are the same as in the general case, the solution of (127) 
takes the form: 



f(r) 



Li(r) 
w 2 W{l) Jo 

Mr) 
w 2 W(l) 

L 3 (r) 

W 2 W(1) J r 

L 4 (r) 



wfW(l) Jo 



dr' r 1 U{r')Wi(r') 

O 

dr' r' U{r')W 2 {r') 

O 

dr' r' U(r')W 3 (r') 
dr' r'U{r')Wi{r'), 



(140) 



and the formulas for the derivatives are analogous to Eq. 
(105)-(107) 

Whenever the VL 2 coefficient vanishes, the recursion re- 
lations established in the previous sections are no longer 



correct. Equations for Ej, i = 1,5, become a system of 
linear equations with a singular coefficients matrix. To 
give an example how to circumvent this problem, let us 
derive an analytical equation for J^£- . In this special case 
Ei and E2 are: 



Ei 



dw 3 



2w{ 



3/(3) 



dfii „„, 



d W3 {r) ~d^ rf W 



' nir a « — / ( r )- fi lfl « — r f\ r ) ( M1 > 

ow 3 ow 3 ow 3 ow 3 

dU(r) 



dw 3 



= 



F 2 df^ df" 1 dSl,. „ 

E 2 = wi -r, 1- fin 1" — / ( r ) 

Ow 3 ow 3 2 ow 3 



1 dQ, 



(142) 



f(r) + V(r)=0. 



2dw 3 

We form the combinations 
d 



E 3 = — (Ei + rEa), 
or 

E 4 = |-(rE 3 -2E 1 ) 
or 



(143) 
(144) 



and solve Eqs. (141)-(144) for |£ instead of -§^. 1 
result is: 

df_ = 

dw 3 



1 

'2^7 

dw 3 
dU"(r) dn 2 



W(r)-2fW + ^r/(r) 
ow 3 Ow 3 



[2/ 



'(r) + r/"(r)]}-^{-8V , (r)-2rV"(r) 



<9fii 



2 ^ + g r/ '' W + S [2/(3)W + r/(4)W] 

(145) 

This result needs now to be formally integrated over r. 
The resulting integrals can be expanded as: 



(146) 



where the superscript (— J) was introduced to denote the 
/-fold integration over r with the boundary condition 
/( _/ )(oo) = (the so-called antidiffcrentiation). The in- 
tegrals of f(r) over r can be obtained by the consecutive 
integration of Eq. (127): 



= 7T <rf"(r) + n x rf(r) - tfH0(r) 



(147) 



The result becomes: 
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df_ 

dw 3 



1 

2£h 



-2M-V(r) + 2V^ 2 \r) - 2 dU [ ^ + ^l[ r f^(r) - / ( " 2) (r)] + ^[/(r) + r/'(r)] 



dw 3 



dw 3 



w 2 



8V(r) - 2rV'(r) + 2V(r) 2^M + §^[r/'(r) - /(r)] + + r/ (3) (r)] 



(148) 



r 



Antiderivatives of the functions U(r) and V(r) can all special case U2 — U3 = W2 — ws = x mentioned earlier 
be obtained in an analytical way. For example, in the they are: 

I 



— U. =_L {-Arw^Uwi + 2x) 2Ei(-rai) - e 4ra Ei(-4rx) - Log(4r) 

OW3 o(wi + 2x) 

+ 2Log(wi + 2x) - 2Log(wi) - Log(x)] + Ae rWl [^rw x {w 1 + 2x) - Ax] + 8(wi + 2x)} , (149) 



V { ~ 1} {r) 



-r{2x+w\ ) 



{-r 2 wie rwi [e irx (2x + wi)Ei(-4rz) - 2(2x + Wl )Ei(-r Wl ) + WiLog(4rx) 

-'fr 2 w 1 (2x + wi) + 4r 2 wi(2x + Wi) tanh -1 ( — - — ) 

\x + Wi J 



2r 2 (2x + w 1 ) 
+ 2xLog(rx) + xlog(16)] + e rwi 
+ Arx + 2] - 2(r(2x + Wl ) + 1)} , 



(150) 



V(- 2 Hr 



iv 1 e 



-r(4x-\-wi) 



( e r{2x+w 1 ) [_ TO 2 e 4rx Ei (_ 4ra .) + 2rw 2 e 2rx Ei[-r (2x + wi)] 



Arx(2x + Wi) 

- 2rxw\e Lrx Ei(— Arx) + Arxwie 2rx Ei[— r(2x + wi)} + rwjLog(4ra;) + ^rw\(2x + wi)] 



+ e 



r(2x-\-wi) 



2rxwiLog(rx) + rxu>iLog(16) — Arw\(2x + wi) tanh 
2rw x {2x + w 1 )Ei(-rw 1 )e r( - 2x+Wl) + Axe 2rx \ . 



( 1 ) 


-Ax 


\ X + Wl J 





r 



(151) 



In a similar way higher-order derivatives over the non- 
linear parameters can be calculated. One needs to use 
the expressions for Ei and E2 differentiated the desired 
number of times over U2, 1*3, 1^2,^3 as a starting point 
and form the same combinations as in the above exam- 
ple. 



B. Vanishing Qi coefficient 

Vanishing fli is a by far less troublesome special case 
than the one considered in the previous subsection. Con- 
ditions under which f2j vanishes are found by recasting 
it into a particular form 



2/2, 2 1 2 1 2\ , / 2 2\l 2 2\ 
h {U 2 + u 3 + W 2 + w 3] + \ U 2 ~ u 3 )( w 3 ~ w 2) 

(152) 



so we may solve fii 
found to be w\ = 



against w\. The result is trivially 



w% ± \fK) , with 



A = (u 2 + u\ + wl + wl) — A(u 2 - u 2 )(w 2 — wl) and 
we see that fii vanishes after some coincidental choice of 
the nonlinear parameters defined by the above equation 
rather than for some particular class of the integrals. In 
this special case the homogeneous differential equation 
takes the form: 

™\rf ( H ] (r) + 2w 2 fH ] (r) + n 2 rf H (r) = (153) 

This equation can be solved by using the recursion rela- 
tions for the coefficients in the series expansions derived 
in the general case, by setting 51i = and recognizing 
the resulting series in terms of the well-known special 
functions. Since we have already presented a detailed 
example of such a procedure, here we only list the fi- 
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nal equations in a convenient form. The solutions of Eq. 
(153) are denoted by double-tilde to distinguish them 
from the previous ones: 



Li(r) = F 3 
L 2 {r) = r F 3 



n 



2 4 

r 



1 3 3 

2' 4' 4' 256^ 

3 5 _Jh_ 

4' '4' 256w? 5 



L 3 (r) =r 2 F 3 



2 4 

r 



5 5 3 

4'4'2' 256w'( 
Jh_ r 
256wl' 





I In 

4 ' 4 ' u 3 



(154) 
(155) 
(156) 
(157) 



where G™" 



ax,..., dp 



,bq 



is the Meijer G-function. 



The solution of the inhomogcneous equation can now for- 
mally be written as: 

f(r) = £l J r) r dr' r' U(r')Wx(r') 



w 2 W(l) Jo 
Mr) 
w 2 W(l) 

wfW(l) Jr 
L 4 (r) 



O 

dr' r'U{r')W 3 {r') 
dr' r' U(r')W 4 (r'), 



(158) 



w 2 W(l) Jo 



with the definitions of W(r) and M^i(r) analogous to Eq. 
(90) and (91), respectively. 

Since the recursion relations derived in the general case 
remain valid for S^i = we can rewrite them as they take 
much simpler form here, for instance: 

+ S [/(r) + rf{r)] ~^ [/ " (r) + r/(3)(r)] } < 

(159) 

so that higher-order derivatives over the nonlinear pa- 
rameters arc calculated from the recursion relation for 
the general case by putting Oi = at the end of each 
recursive step. 



C. Vanishing Qj and Q 2 coefficients 

Situation when Hi = and f2 2 = is quite rare since 
the conditions given in the two previous subsections that 
make f2i and f2 2 vanish must mutually be satisfied. This 
occurs, for example, when u 2 = u 3 = u> 2 = w 3 = x 
and additionally Wx — 2x. The homogeneous differential 



equation has disarmingly simple four linearly indepen- 
dent solutions: 



Lx(r) = 1, 

L 2 {r) = r, 

L 3 (r) = r 2 , 

Li(r) = r Log(r) — r. 



(160) 
(161) 
(162) 
(163) 



The above solutions were denoted by check mark to sepa- 
rate them from the previous ones. The Wronskian W(r) 
and the Wi(r) determinants can be brought into the fol- 
lowing closed forms: 

(164) 
(165) 

)g(r), (166) 
(167) 
(168) 

so that the solution of the inhomogeneous differential 
equation takes the form: 



W{r) = 


2 


Wx(r) = 


r, 


W 2 {r) = 


2 - 




1 


W 3 (r) = 




r 


Wx{r) = 


-2, 



f(r) 



1 



UK 



dr' r' 2 U(r') + r dr' r'U(r') [1 - Log(r)] 



dr' U{r') + r [Log(r) - 1] / dr' r'U{r') 

Jo 

(169) 

where for example in the case u 2 = u 3 = w 2 = w 3 = x, 
Wx = 2x: 



U(r) 



4xe 



{2r 2 x 2 e 2rx [e 4ra Ei(-4ra-) - 2Ei(-2ra;) 



+ Log(ra;)] + 6rx + e 2rx [2rx(jrx - 1) - 1] + 1} 



(170) 



However, even in such a simple case not all of the above 
integrals can be calculated fully analytically, so we still 
need to struggle with the numerical integration. A lit- 
tle bit more difficult is the differentiation of the master 
integral with respect to the nonlinear parameters. For 
example, the derivative over w 3 is obtained from the spe- 
cial forms of the two identities which were derived in the 
previous subsections: 



El — — LU 1 

<9Qi 



Wi r — 2w 1 — (r) 



dw 3 dw 3 



r/"(r)-fV(r) 
aw 3 ow 3 

a^2 ( , | dU(r) =Q 

dw 3 dw 3 



E 2= ^a/^ + iafii 

dw 3 2 dw 3 



(171) 



(172) 
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We take the combination (Ei + rE^) to cancel out the 
term 



df^ ill i • • • df^ 

q Ws and solve the resulting equation against g 



1 



dU{r) 
dw* 



'2d^ rf (r) -2a^ r/(r) 



This equation needs now to be antidifferentiated three 
times to give: 



df_ 

dw 3 



1 

dSli 
dw 3 



5C/(- 3 )(r) 



dw 3 



5 ^[r/t" 8) (r) -/<"-» (r)] 



(174) 



The antiderivatives of the master integral are obtained 
by consecutive antidifferentiation of the expression: 



w 2 rf^(r) + 2w\fW (r) = U(r), 



e -9- f 



rf(r) 



l/(-4)( r ) 



(175) 



Higher-order derivatives with respect to the nonlinear pa- 
rameters are obtained using the same procedure albeit Ei 
and E2 need to be differentiated an arbitrary number of 
times with respect to 1*2,^3, W2,w 3 before putting Oi = 
and £^2 = 0. 



D. Vanishing f)f — 4tof Q2 coefficient 

In the special case Cl 2 = two pairs of roots 

of the a polynomial, Eqs. (25) and (31), lying on 
the same side of the complex plane coincide, so that 

(7 = w\ It 2 + 2?) ■ As a result, the homogeneous dif- 



ferential equation, Eq. (54), can be brought into the 
form: 

w 2 A [ri/ H (r)] = 0, (176) 
where A is a differential operator defined as: 

i= |5>-^ ( i77 ) 

with q 2 = —^5-. Eq. (176) can be solved by decom- 
posing it into a system of two second-order differential 
equations: 



w\A [rh(r)} = 0, 
Af H (r) = h(r). 



(178) 
(179) 



The first of these equations has the form rh"(r) + 2h'(r)- 
rq 2 h(r) = 0, so that the general solution is: 



-■<!'■ 



-qr 



h( r )=C 1 —+C 2 



(173) an d Eq. (179) takes the form: 



fUr)-1 I f H {r)=C l — 
r 



C 2 



-qr 



(180) 



(181) 



The latter equation is solved with elementary methods. 
Finally we conclude that Eq. (176) has four linearly in- 
dependent solutions that can be chosen as: 



e qr , e- qr Ei[2qr] - e qr Log[2qr], 



and e qr Ei[-2qr] 



-qr 



Log[2gr]. 



(182) 



The solution of the inhomogeneous differential equation 
takes the form analogous to Eq. (101). Similarly, Eqs. 
(105)— (107) are the derivatives of the master integral 
with respect to r. 

Since Dq = Q, 2 — 4w 2 fl2 appears in the denominator in 
nearly all recursion relations derived for the general case 
they become invalid here. However, this problem can be 
circumvented by using the same trick as in the fl 2 = 
case, namely solving the system of Eqs. (141)— (144) with 
respect to Jj^- and performing consecutive antidifferenti- 
ations. Since the derivation is exactly the same as in the 
subsection VII A there is little point in repeating it here. 



VIII. NUMERICAL EXAMPLES 

In this section we present results of calculations on the 
representative set of master integrals with some hand- 
picked values of the nonlinear parameters. We imple- 
mented a general code that is able to calculate the values 
of the master integral with arbitrarily chosen nonlinear 
parameters. The code is written in the C programming 
language and all the calculations were performed in the 
quadruple arithmetic precision using the GCC Libquad- 
math library. Handful of the results presented here were 
additionally checked by using an independent program 
written in Mathematica with the octuple arithmetic pre- 
cision. Comparison with the results obtained in the ex- 
tended precision shows that calculations performed in 
quadruple-precision, using 101 points of the Tanh-Sinh 
quadrature [126, 127] for all numerical integrations, gave 
an accuracy of at least long double precision (around 20 
significant digits) and much better on the average. 

It is easy to verify that for the calculation of the mas- 
ter integral one can also use a different procedure, based 
on the series expansion of exp(— w±r 12) around Wi = 
under the sign of the integral in Eq. (8). Since the latter 
expansion is uniformly convergent for any positive value 
of wi one can perform term by term integration what 
leads to the identity: 



f (r; Wl ) = J2fn(r;0) 
n=0 



(-l) n w[ 



(183) 
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TABLE I: Comparison of the values of the master integral 
with U2 — W2 = 1, U3 = W3 = 2, wi = 1.0 calculated ac- 
cording to Eqs. (176) and (101). Calculations performed for 
r = 2. [k] denotes 10 fe . 



expansion length 


f(r) value 




10 


2.38528323813123779081[- 


-04] 


20 


2.39121157642870208323[- 


-04] 


30 


2.39121170158061051314[- 


-04] 


40 


2.39121170158262980140[- 


-04] 


50 


2.39121170158262983284[- 


-04] 


60 


2.39121170158262983284[- 


-04] 


70 


2.39121170158262983284[- 


-04] 


75 


2.39121170158262983284[- 


-04] 


Eq. (101) 


2.39121170158262983284[- 


-04] 



TABLE III: Comparison of the values of the master integral 
with U2 = W2 = 1, U3 = W3 = 2, w\ = 2.0 calculated ac- 
cording to Eqs. (176) and (101). Calculations performed for 
r = 2. [k] denotes 10 fc . 



expansion length 


f(r) value 




10 


-3.63011904625504245573[- 


-04] 


20 


1.07042192138901654774[- 


-04] 


30 


1.17367132416023521971[- 


-04] 


40 


1.17538021409679753524[- 


-04] 


50 


1 . 17540746158661903854[- 


-04] 


60 


1 .17540789464360697759 [- 


-04] 


70 


1.17540790155566410929[- 


-04] 


75 


1 .17540790165897003951 [- 


-04] 


Eq. (101) 


1 . 17540790166845070630[- 


-04] 



TABLE II: Comparison of the values of the master integral 
with U2 — W2 = 1, it3 = 1U3 —2, wi = 1.5 calculated ac- 
cording to Eqs. (176) and (101). Calculations performed for 
r = 2. [k] denotes 10 fe . 



expansion length 


f(r) value 




10 


1.32952081604592501191[- 


-04] 


20 


1.63128202122143112608[- 


-04] 


30 


1.63165160669687214817[- 


-04] 


40 


1.63165195079171515995[- 


-04] 


50 


1.63165195110063806339[- 


-04] 


60 


1.63165195110091457699[- 


-04] 


70 


1.63165195110091482555[- 


-04] 


75 


1.63165195110091482578[- 


-04] 


Eq. (101) 


1.63165195110091482597]- 


-04] 



TABLE IV: Comparison of the values of the master integral 
with U2 = W2 = 1, U3 = ins = 2, wi = 2.5 calculated ac- 
cording to Eqs. (176) and (101). Calculations performed for 
r = 2. [k] denotes 10 fe . 



expansion length 


f(r) value 




10 


-3.96382454835359866193[ 


-03] 


20 


-7.39453380391436518395[ 


-04] 


30 


-3.92253435711155687731[ 


-05] 


40 


6 . 940930 190407504 1 6082 [ 


-05] 


50 


8.55501340987912083454[ 


-05] 


60 


8.79376564476664924936[ 


-05] 


70 


8.82922222853990351131[ 


-05] 


75 


8.83787155222436438480[ 


-05] 


Eq. (101) 


8.83545586039834446778[ 


-05] 



The above series is convergent for any value of w\ and 
gives exactly the same numerical result as Eq. (101). 
The prescription how to calculate the integrals with an 
arbitrary power of r%2 but wi = 0, f n {r\ 0), was recently 
presented by means of the open-ended recursion relation, 
c.f. Eq. (48) of Ref. [44]. Therefore, Eq. (176) is an in- 
teresting alternative to the analytical equation derived in 
the previous subsection. It is worth considering in details 
how fast the above series expansion converges for a given 
value of wi and how many terms are necessary to obtain 
long double precision which one can easily get by using 
Eq. (101) throughout. To make such a comparison pos- 
sible, we implemented the mentioned recursion relation 
to advance the power of r±2 in / n (^;0) as much as nec- 
essary in the Mathematica package. However, the first 
problem encountered was the numerical stability of this 
recursion. Although the starting values for the recursion 
were computed in the octuple arithmetic precision, after 
n — 75 steps only few digits were estimated to be correct. 
Therefore, it is rather pointless to go beyond this value 
of n. On the other hand, further extension of the arith- 
metic precision which is already two times bigger than for 
calculations based on Eq. (101) will slow down the code 
dramatically and make it inferior to the numerical inte- 
gration approach. In Tables 1-5 we present a comparison 
of the values of the master integral obtained according 



to the Eq. (176) with different expansion lengths and 
obtained by numerical integration in Eq. (101). We have 
chosen representative values of the nonlinear parameters 
(v,2 = W2 = I, Us = W3 = 2) which were kept fixed and 
we have varied the value of w\ to examine the behaviour 
of the expansion Eq. (176) with increasing w\. We see 
that for W\ = 1 series expansion defined by Eq. (176) 
converges fast and smoothly towards the correct value 
and only a few tens of terms are necessary to obtain the 



TABLE V: Comparison of the values of the master integral 
with U2 = W2 = 1, M3 = W3 = 2, w\ = 3.0 calculated ac- 
cording to Eqs. (176) and (101). Calculations performed for 
r = 2. [k] denotes 10 fc . 



expansion length 


f(r) value 




10 


-2.28509213221681747900[ 


-02] 


20 


-2.90220438815452281770[ 


-02] 


30 


-2.77016299506271493255[ 


-02] 


40 


-2.54662563339713840590[ 


-02] 


50 


-2.33351029802116868289[ 


-02] 


60 


-2.14729270210488323595[ 


-02] 


70 


-2.19874109964523467927[ 


-02] 


75 


-2.19299042315798720906[ 


-02] 


Eq. (101) 


3.89692438286734239005[ 


-06] 
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long double precision result. An even better behaviour is 
met for lower values of W\. However, when W\ increases 
beyond one the convergence of Eq. (176) deteriorates and 
even for w\ = 1.5 as much as 75 terms of the expansion 
are not enough to obtain a reliable precision of 21 sig- 
nificant digits. For w\ = 2.5 only three significant digits 
are recovered after 75 terms and for wi = 3.0 the series 
converges so badly that no useful information about the 
value of the master integral is obtained after 75 terms. 
Moreover, the result for w\ = 3.0 is clearly wrong since 
by a simple inspection of Eq. (8) we observe that the 
value of the master integral is always positive. 

Although a simple comparison provided in the above 
clearly shows that the numerical integration approach is 
superior compared to the series expansion method, we 
must admit that the numerical integration has its own 
problems. They appear for small values of wi, say, lower 
than 0.2. In this regime, the integrands in Eq. (101) 
vanish slowly and significant contribution to the value 
of the master integral comes from the large r' in inte- 
gration over the interval [r, +oo [. The treatment of such 
situations requires an efficient matching of the series ex- 
pansion around with the asymptotic expansions of the 
functions Lj(r). Moreover, for large r' accurate calcula- 
tion of Wi(r') becomes difficult because significant loss 
of digits occurs due to the subtraction of two near-equal 
numbers. Therefore, we believe that the most efficient 
method of calculation of the master integral will be a 
suitable union of two algorithms described here. Series 
expansion is to be used for small w± where it converges 
fast and only a handful of terms is required to obtain 
desired accuracy. For larger w\ numerical integration is 
superior and is able to provide arbitrary accuracy. In Ta- 
ble 6 we additionally listed values of the master integral 
with some combinations of the nonlinear parameters cor- 
responding to the general case without comparing them 
to the series expansion method. 

Special cases of the master integral were implemented 
separately taking advantage of the fact that Li(r) func- 
tions are expressed in terms of known special functions. 
All necessary special functions were implemented using 
Chcbyshev interpolation method. In Table 7 we give ex- 
amples of the values of the master integral in one impor- 
tant special case corresponding to the exponential ver- 
sion of the symmetric James- Coolidge basis set, namely 
U2 = U3 = W2 = W3 = x with an arbitrary value of w\ 
and x (vanishing coefficient). 



IX. OUTLINE FOR THE FUTURE 

In this paper we introduced a new explicitly correlated 
basis set for state-of-the-art ab-initio calculations on di- 
atomic molecules and reported analytical formulas ready 
to apply for all molecular integrals appearing in the non- 
rclativistic calculations. While a physical application of 
the theory presented here will be reported soon, we would 
like to stress that our theoretical results will find several 



TABLE VI: Examples of the values of the master integral 
with U2 = W2 = x, us — «)3 = y calculated according to Eq. 
(101). The symbol [k] denotes the powers of 10, 10 fc . 





X 


a 
y 


T 


/(r) value 




2.0 


2.0 


3.0 


2.0 


1.20162929654132118557[- 


06] 


2.0 


2.0 


3.0 


6.0 


1.85690183473424018070[- 


-15] 


3.0 


2.0 


3.0 


4.0 


2.54651823562870024818[- 


■11] 


3.0 


2.0 


3.0 


1.0 


1.27020704100393024656[- 


■04] 


2.0 


2.5 


1.5 


5.0 


5.46845126485930791142[- 


-11] 


2.0 


2.5 


1.5 


10.0 


1.23085262413134222873[- 


■11] 


5.0 


1.0 


3.5 


1.0 


1.34862879254745038343[- 


-04] 


5.0 


1.0 


3.5 


3.0 


1.12445814382106648815[- 


-08] 


8.0 


2.0 


4.0 


1.0 


1.12012736029631067277[- 


-05] 


8.0 


2.0 


4.0 


0.5 


2.73482762158894864378[- 


-04] 



TABLE VII: Examples of the values of the master integral 
with U2 — ui2 = «3 = W3 = x, calculated according to Eq. 
(140). The symbol [k] denotes the powers of 10, 10 fc . 





X 


r 


/(r) value 




1.0 


1.0 


1.0 


1.7031252809284U22796[ 


-02] 


1.0 


2.0 


1.0 


1.014022U417374426756[ 


-03] 


1.0 


3.0 


1.0 


8.004150U434536226469[ 


-05] 


3.0 


1.0 


1.0 


5.95180137043458120693[ 


-03] 


3.0 


2.0 


1.0 


4.27403013549289339748[ 


-04] 


3.0 


3.0 


1.0 


3.72976454960392088836[ 


-05] 


1.0 


1.0 


6.0 


7.05646817317973364415[ 


-07] 


1.0 


2.0 


6.0 


1.84467129942180429U8[ 


-12] 


1.0 


3.0 


6.0 


6.55011946233785000932[ 


-18] 


3.0 


1.0 


6.0 


1.43442768276386167206[ 


-07] 


3.0 


2.0 


6.0 


4.47644083U8078530669[ 


-13] 


3.0 


3.0 


6.0 


1.76871016217927283912[ 


-18] 



important applications. 

First of all, the Slater gcminal basis is expected to im- 
prove the convergence of molecular calculations on two- 
electron diatomic molecules by several orders of magni- 
tude. With the advent of new experimental tools that 
allow measurements of the dissociation energy of H2 with 
an astonishing accuracy of 10~ 4 cm -1 [128] and with the 
announcements that this level of accuracy will be im- 
proved by two orders of magnitude, new molecular cal- 
culations will be necessary to reproduce the experimen- 
tal data. Especially important in this respect will be the 
calculation of the relativistic integrals in the basis of the 
Slater geminals. We expect that the accuracy of the rel- 
ativistic corrections reported in Rcf. [129], computed in 
the basis of explicitly correlated Gaussian geminals, can 
be improved by a few orders of magnitude. Also the QED 
effects could be accounted for in a more accurate way, to 
produce not only state-of-the-art estimates of the disso- 
ciation energy, but also of the rotational and vibrational 
spacings [130]. 

The second important application of the theory pre- 
sented above is the numerical calculation of the inte- 
grals in the basis set of Slater orbitals for diatomic 
molecules. While the theoretical background was intro- 
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duced by Pachucki in 2009 [44] , his algorithms for certain 
classes of integrals turned out to be inefficient for prac- 
tical implementation. Using the geminal recursion rela- 
tions and putting the exponent in the correlation factor 
Wi equal to zero, one obtains much simpler and numeri- 
cally more convenient recursion relations. 

Also worth mentioning are the calculations of the rel- 
ativistic integrals in the basis set of the Slater orbitals 
for diatomic molecules. At present, no ab initio program 
for molecular calculations has available all integrals ap- 
pearing in the Breit-Pauli theory, even in the Gaussian 
basis set. We plan to apply our theory to the expres- 
sions for the most difficult class involving the rf 2 2 factor, 
and perform actual calculations with just one numerical 
integration in one dimension. In this way, accurate cal- 
culations of the fine and hypcrfine structure of diatomics 
will become possible. 

Finally, the basis set of the Slater geminals can be 



used in the explicitly correlated MBPT/CC theories, 
thus greatly improving the accuracy of the present ap- 
proaches based on the Gaussian orbitals and linear or 
exponential correlation factors. 
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Appendix A: Matrix elements of the kinetic energy 
operator 

Let us consider the two-center two-electron Schodinger 
Hamiltonian: 



1. 



1, 



Zr 



Zr 



1 



1 



H = --V! - -Vi - — - — - — - — + — + - 

2 2 via r 2 A t\b r 2B r 12 r 

(Al) 

where Zk denotes the nuclear charge of the nucleus K 
and the notation for the other quantities is the same as 
in Eq. (7). The basis functions are of the form (6). The 
overlap integrals between these basis functions and the 
matrix elements of the nuclear attraction and electronic 
repulsion operators are obviously expressed through the 
integrals from the family (7) with proper powers of TiK 
and ri2- The only difficulty is to express the matrix ele- 
ments of the kinetic energy operator through the integrals 
(7). Let us introduce a shorthand notation that will be 
used throughout this Appendix: 



Our derivation was inspired by the procedure given by 
Kolos et al. [25, 26]. In fact their result is the w\ = 
w[ = limit of our equation. The matrix element is 
transformed as: 



(ijkln\ - V x |i j kin) = J d r 2 <£fc;( 2 ) tpk'i>(2) 

(A3) 

where the Green's theorem was used. Let us now consider 
only the integration over the coordinates of the first elec- 
tron. For simplicity wc will consider the case n = n' = 0. 
The general form of the matrix element will then be ob- 
tained by differentiating n times with respect to wi and 
n' times with respect to w[ , and multiplying by the factor 
(— 1)"+" . The latter integral is rewritten as: 



d 3 r, V x [^(l)e-™] - Vi [^(lK^ 
d 3 n { [Vi^ (1)] [Vi^v (i)] e -(«'i-K)r ia 

T K[V!^(i)]^(i) 



(A4) 



+^i^(l)[V 1 ^v(l)])[Vie-^+^)'--]} 
d 3 n [ Wl w' mj (l) Wf (l)e-^ +w >^ 



-tW[Aiw(1)]w(1) 

+w lifij (l)[A l <p i , f (l)})e-^ +w >"}, 

where the Green's theorem was used in the last step to 
remove the gradient operator working on terms contain- 
ing the ri2 factor. The equation for the Laplacian of the 
Vij(l) is rather straightforward to derive and the final 
result is: 



A mj {l) = i(i + j + 1)^_ 2)J -(1) + j(i + j + % w _ 2 (l) 



+ («§ + u|) ip^l) - u 3 ( i + | + 2 ) ^-1^(1) 



- "2 ( j + - + 2j tfi,j-i(l) - r ijtpi- 2 ,j-2(l) 

+ -iu 2 rVi-2,j-i(l) + 2JU3r 2 <Pi-i,j-2(l) 

+ u 2 u 3 [<^ + ij_i(l) + <^_i, i+ i(l) - r 2 y>i_i i; ,-_i(l)] 



1 . , 1 

-lU 2 <Pi-2,j+l{l) - -J« 3 ^+l,3-2(l) 



\ijkln) 



,.j 



j 



' 1A' \B'1A' 2B' 12 c 
yii(l)Va( 2 )^12 e ~ 



-U3riA—U2riB —W2r2A — w 3 r2B—w 1 r 1 2 



(A5) 



(A2) 



By inserting the formulas (A4) and (A5) into Eq. (A3) 
one arrives at: 
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(ijklO\ - V 2 \i'j'k'l'0) = Wiw'i(ijklO\i'j'k'l'0) 



-^-T W + J +!)(«- 2, jklO\i'j'k'l'0) + j(i + j + - 2, klO\i'j'k'l'0) 

+ (ul + u\) (ijklO\i'j'k'l'0) - u 3 (i + 3 - + 2^j (i - l,jklO\i'j'k'l'0) - u 2 (j + | + 2^j (ij - 1, klO\i'j'k'l'0) 

~r 2 ij(i~2,j - 2,kl0\i'j'k'l'0) + ^iu 2 r 2 (i-2,j - l,klO\i'j'k'l'0) + ^ju 3 r 2 (i- 1, j - 2, klO\i'j' k'l'O) 
+u 2 u 3 (i + l,j - l,klO\i'j'k'l'0) + u 2 u 3 (i - l,j + l,MO\i'j' k'l'O) - r 2 u 2 u 3 (i - l,j - l,klO\i'j'k'l'0) 

--iu 2 (i - 2, j + 1, klO\i'j'k'l'0) - ^ju 3 (i + 1, j - 2, klO\i'j'k'l'0) 



(A6) 



Wl 



wi + w\ 



T [*'(*' + f + l)(ijklQ\i' - 2,j'k'l'0) + MO\i'j' - 2, k'l'O) 



+ K 2 + (ijklO\i'j' k'l'O) - u' 3 ( i' + - + 2 ) (ijklO\i' - I, j 1 k'l'O) - u' 2 [ j' + - + 2 ) {ijklO\i'f - 1, k'l'O) 



r 2 i'j'(ijklO\i' -2,j' - 2, k'l'O) + -i'u' 2 r 2 (ijklO\i' -2,f - 1, k'l'O) + -j'u' 3 r 2 (ijklO\i' - l,j' - 2, k'l'O) 



1 

2" ~ z ' ' ' 2 

+u' 2 U3(ijfcZ0|i' + 1,/ - l,fc'/'0) +u' 2 u' 3 (ijklO\i' - 1, / + l,fc'Z'0) - r 2 i4?4(ijfcZ0|i' - 1, j' - 1,/c'Z'O) 

-iiV^ijWOli' - 2, j' + 1, k'l'O) - \j'u' 3 {ijklO\i' + l,f - 2, k'l'O) 



where the notation (ijklO\i'j' k'l'O) was used to designate 
ordinary overlap integrals which belong to the integral 
family (7). The above equation is already symmetric with 
respect to the interchange of primed and non-primed in- 
dices. It is noteworthy that the singularity appearing 
when wi = w'i = is only apparent. It can easily be re- 
moved by taking first the limit w[ — > W\ and then putting 
W\ =0. To obtain the most general matrix element of the 
kinetic energy operator, (ijkln\ — S7 2 \i'j'k'l'n'), one needs 
to differentiate Eq. (A6) with respect to w\ and w[. This 
differentiation is easily carried out explicitly with any 
symbolic mathematical program. The resulting formulas 
are quite compact and make their direct implementation 
straightforward. On the other hand, by multiplying both 
sides of Eq. (A6) by wi + w[ and performing consecutive 
differentiation, it is quite easy to derive a recursion re- 
lation that connects the values of (ijkln\ — Vf |i' 'j'k'l'n') 
with different n and n' . However, this recursion is inher- 
ently unstable with the increasing n and n', and can be 
considered inferior to the approach based on the analyt- 
ical equations, at least for larger n and n' . 



Appendix B: Analytical formulas for the functions 

W(r) and V(r) 



In this Appendix we list explicit analytical for- 
mulas for the functions W(r; w±, u 2 , u 3 , w 2 , w 3 ) and 
V(r; Wi, u 2 , u 3 , w 2 , w 3 ) which are defined as the inverse 
Laplace transforms in Eqs. (52) and (53), respectively. 
We would like to stress that in the formulas given below 
all the terms proportional to the Dirac delta distribution 
or its derivatives were omitted. This can be done because 
for r > these terms never contribute to the values of 
the master integral derivatives calculated from the recur- 
sion relations presented in the text. However, the missing 
terms can be recovered by taking the Laplace transform 
of the listed equations and comparing with the proper 
analogue of Eq. (26). This can easily be done by using 
any symbolic mathematical package. 

The function W(r) takes the form: 



6 

W(r) = Y,Wi(r) 

i=l 



(Bl) 
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with 



W 2 (r) = Wl e- r( - U3 + w i+ w ^ LL + 



1 U2 + W1+W2 



1 U 3 + Wi + w 3 



(B2) 
(B3) 



W 3 (r) 



WlW2{v.2 — w 3 ) + U 2 W 2 ~ U 3 W 3 1 

(lUl + W2) 2 — w 3 r 
Wi(wi + w 2 ) 1 



|_ e -r(« 2 +«, 3 ) [ 2 + 2r(u 2 + ttf 3 ) + r 2 (u 2 + w 3 f 



+ e - r ^ +W2+W1 ^ [2 + 2r(u 2 + w 2 + w{) + r 2 (u 2 + w 2 + wi) 2 ] } , 



(B4) 



Wi{r) = 



WiW 3 (u 3 — W2) + UgWg — 1 



(wi + w 3 ) 2 — w\ 



D -r(u 3 +w 3 +wi) 



-r(u 3 +W2) 



Wl(wi + w 3 ) 1 
(wi + w 3 ) 2 — w\ r 3 



|_g-r(« 3+W2 ) [ 2 + 2r ( U3 + W2 ) + r 2 (u3 + W2 y 



-r(u 3 +w 3 +w 



l} [2 + 2r(u 3 +w 3 + Wl ) + r 2 (u 3 + w 3 + Wl f] X 



(B5) 



W 5 (r) 



wiu 2 {wl - m|) + - ugwg 1 
(it 2 + wi) 2 - u\ r 



-r(v,2+W2+Wl) 



-r(u 3 +w 2 ) 



jg-r(« 3 +^) [ 2 + 2r ( W3 + W2) + r 2( M3 + 



11% — U 2 + U2W1 1 

(«2 + wi) 2 — u 2 r 3 

_ e -r(u 2+W2+Wl ) [ 2 + + W2 + Wl ) + r 2 (M2 +W2+ Wi) 2] J 



(B6) 



w w = wiu 3 {u 2 - w 3 ) + u 2 w 2 - u 3 w 3 1 f_ r{u3+W3+Wl) _ e - r{u2+W3) 
-(u 2 + wi) 2 + u\ r I 



Uo — vM — U 3 Wi 1 



|g-r(« 3+IU3 ) [ 2 + 2r(u 2 + life) + r 2 {u 2 + w 3 f 



— (u 3 + Wi) 2 + u\ r 3 

_ e -r(u 3+W3+Wl ) [ 2 + 2r(us +W3+ + r 2 (u3 + ^ + Wi) 2] J 

I 



(B7) 



Let us denote the permutation U2 -H- W3, U3 O 1^2 by 
P12 and the permutation U2 O U3, w 2 O W3 by Pab- 
The reason for adopting such a notation becomes clear 
when one considers the symmetries of the master inte- 
gral. Calculation of W(r) is simplified by the following 
relations: 

P AB W 1 (r) = W 2 (r), 

P AB W 3 {r) = Wi(r), (B8) 

PA B W 5 {r) = W 6 (r), 

so that the programming effort is halved. Further sim- 
plifications occur after observing, for instance, that: 

e -r(u 3 +w 2 ) 

- [2 + 2r(u 3 + w 2 ) + r 2 (u 3 + w 2 f] = 

Q2 [ p -r{u 3 +w 2 )l ( B9 ) 

dr 2 



so that in the implementation one can concentrate on 
the calculation of the derivatives of the quantities like 
e -ar j r smce an y derivative of W(r) with respect to the 
nonlinear parameters and r is expressed through them. 

Explicit form of V(r) is expressed conveniently as: 



V{r) = Vx{r) + V 2 {r) + V 3 (r) + ciV 4 (r) + c 2 V 5 (r), 

(BIO) 
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with 



and 



e — r(u 3 +w 2 +w 1 ) „— t(v.3+W2) 

Vi(r) = (w\ - u% - u 3 w 3 ) H (-10? + U3 - 



u 3 w 3 - 



r 

,2 „„2\ 



3 — r(jI 2 +U>2+M>i) 



U2W3- 



V a (r) = 



WiW 3 (u 3 — W$) + — U2W2 1 

(wi + w 3 ) 2 — u>2 r 
w\ (wi + w 3 ) 1 



g-r{«s+««3+Wi) _ g— r(«3+tU2) 



|_g-r(« 3 +^2) [ 2 + 2r (u 3 + w 2 ) + r 2 (u 3 + w 2 f 



(wi + w 3 ) 2 — w\ r 3 
+ e - r ("3+-3+»i) [ 2 + 2r(u 3 + w 3 + wi) + r 2 (u 3 + w 3 + wif] } 



V 3 (r) 



Wiiu 3 (w 2 - U§) + W2W 3 (u 3 — u\ + w\) 1 



(wi + w 2 ) 2 - w 3 



e -r(u 2 +W2+w 1 ) _ g— r(u3+wa) 



WlW 3 



\ {-e~ r ( U2+w ^ [2 + 2r{u 2 + w 3 ) + r 2 (u 2 + w 3 f 



(wi + w 2 ) 2 — w 3 r 3 
e -r(u 2+W2+Wl ) ^ + 2r (u 2 + w 2 + Wl ) + r 2 {u 2 + w 2 + Wl ) 2 ] } , 

V 4 (r) = e -K«2+™ 3 ) {-Ei [-r (wi + u 3 - u 2 )] + Ei [-r (u 3 - u 2 - w 3 + w 2 ) 
- Ei [-r ( Wl - w 3 + w 2 )}} - e r{u2+Wi) B\ [-r (u 2 + u 3 + w 2 + w 3 )} 



+e -r(u 2 +w 3 ) L 



(w\ +w 2 - w 3 ) (wi - u 2 + u 3 ) (u 2 + u 3 + w 2 + w 3 ) 



(wi +w 2 + w 3 ) (u 2 + u 3 + wi) (u 2 - u 3 - w 2 + w 3 ) 
V 5 (r) = e -(-3-«2) Ei [_ r (wi + u , 2 + W3 )j + e r(u 2 - W3 ) Ei [_ r (wi +U2 + u ^ 



Cl = - [u 2 (iwf - W2 + W3) + ^3 («2 - u l + w l)] 
C 2 = i [«2 (w? - W2 + «>3) - W 3 («2 - u 3 + u 'l)] 



(Bll) 



(B12) 



(B13) 



(B14) 

(B15) 

(B16) 
(B17) 



We see that the form of V 2 (r) and V 3 (r) is analogous 
to the Wi(r), i = 3,6. Similarly, 14(r) and Vs(r) are 
expressed through the same combinations of functions as 
U 3 (r) and U 2 (r), respectively. The only difference is that 
some terms contribute with the sign reversed. Since the 
form of Ui(r) is relatively simple and straightforward to 
implement, arbitrary derivatives of the function V(r) can 
be computed by using a proper union of the algorithms 
for U(r) and W(r) functions. 



Appendix C: Auxiliary recursion relations 

In this Appendix we list formulas for the derivatives 
of and jr^- over r up to the third-order which re- 

0W3 awi 1 

suit from the solution of the set of equations for Ej and 
Ei, i = 1,5. Higher-order derivatives can be obtained 
by successive differentiation of the geminal differential 
equation: 



df_ 

dw 3 



2 {Aw 2 n 2 - ilf 

2 

w\ 

+ 4w 2 n 2 - n 2 



dili, 



-2rV{r) 
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dw 3 
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(4)/ 



dw 3 
dU"(r 
dw 3 
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dw 3 
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